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ABSTRACT 

A  computer  program,  called  the  Planetary  Ephemerls  Program  (PEP),  is  being  written 
at  Lincoln  Laboratory.  The  purpose  of  the  program  is  to  improve  planetary  and  lunar 
ephemerides  using  the  results  of  radar  and  optical  observations.  In  this  report,  we 
derive  the  differential  equations  that  are  numerically  integrated  in  PEP  to  determine 
as  functions  of  time  the  positions  and  velocities  of  the  planets,  of  the  Earth-Moon 
barycenter  and  of  the  Moon,  and  the  partial  derivatives  ot  these  positions  and  ve¬ 
locities  with  respect  to  initial  conditions,  masses  and  other  parameters.  Newtonian 
theory  with  the  usual  unrigorous  general  relativistic  corrections  is  employed.  The 
equations  of  motion  and  the  equations  for  the  partial  derivatives  with  respect  to  ini¬ 
tial  conditions  are  presented  in  the  form  needed  in  the  Encke's  method  cr  integration 
used  in  PEP. 
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GENERATION  OF  PLANETARY  EPHEMERIDES 
ON  AN  ELECTRONIC  COMPUTER 


I.,  INTRODUCTION 

A  computer  program,  called  the  Planetary  Ephemeris  Program  (PEP),  is  being  written  at 
Lincoln  Laboratory.  The  purpose  of  the  program  is  to  improve  planetary  and  lunar  ephemerides 
using  the  results  of  radar  and  optical  observations.  The  procedure  for  improving  ephemerides 
is  as  follows.  First,  we  integrate  the  differential  equations  of  motion  of  ihe  planets,  of  the  Earth- 
Moon  barycenter  and  of  the  Moon  using  provisional  values  for  the  various  parameters  (such  as 
initial  conditions  and  masses)  appearing  in  the  theory  of  gravitation  and  motion  employed.  We 
also  integrate  the  differential  equations  for  the  partial  derivatives  of  the  position  and  velocity  of 
the  bodies  with  respect  to  these  parameters.  Then,  for  each  radar  and  optical  observation,  we 
calculate  the  theoretical  values  of  the  measurements  made  in  the  observation  and  the  partial 
derivatives  of  these  theoretical  values  with  respect  to  the  parameters.  Using  the  differences  be¬ 
tween  the  observed  and  theoretical  values  of  the  measurements,  the  errors  of  the  measurements 
and  the  partial  derivatives  of  the  theoretical  values  of  the  measurements,  we  form  the  normal 
equations  and  solve  them  to  get  corrections  to  the  parameters.  With  the  adjusted  parameters, 
we  reintegrate  the  equations  of  motion  and  the  equations  for  the  partial  derivatives;  applying 
the  results  of  these  integrations,  we  again  form  the  normal  equations  and  solve  them  to  get  fur¬ 
ther  corrections  to  the  parameters.  This  process  is  repeated  until  we  obtain  convergence.  Us¬ 
ing  the  parameters  thus  obtained  in  the  integration  of  the  equations  of  motion,  we  generate  ephe¬ 
merides  which  best  agree  with  the  observations  in  a  least-squares  sense. 

In  this  report,  we  present  derivations  of  the  foi  mulas  used  in  PEP  to  determine  as  functions 
of  time  the  positions  and  velocities  of  the  planets,  of  the  Earth-Moon  barycenter  and  of  the  Moon, 
and  the  partial  derivatives  of  these  positions  and  velocities  with  respect  to  the  various  param¬ 
eters.  In  a  second  report,  we  will  describe  the  derivations  of  the  formulas  used  in  the  compari¬ 
son  of  theory  and  observation  and  in  the  least-squares  analysis  parts  of  the  program.  In  a  third 
report,  we  will  present  the  documentation  of  the  computer  program  itself. 

In  deriving  the  equations  in  this  report,  we  employ  Newton's  theory  of  gravitation  and  motion 
with  the  usual  unrigorous  general  relativistic  corrections.  Since  one  of  our  purposes  in  analyz¬ 
ing  radar  and  optical  observations  with  PEP  is  to  test  general  relativity,  the  equations  of  motion 
employed  should  be  derived  in  strict  accordance  with  a.  principles  of  this  theory.  As  explained 
in  Sec.iV-C.  the  equations  of  motion  of  the  general  relativistic  N-body  problem  have  been  derived 
in  principle;  it  is  only  necessary  to  learn  the  derivation  and  to  put  the  equations  in  a  form  amen¬ 
able  to  the  generation  of  ephemerides.  The  equations  will  resemble  the  Newtonian  equations 
with  rigorous  rather  than  unrigorous  general  relativistic  corrections.  Tne  rigorous  general 
relativistic  treatment  will  be  presented  in  a  subsequent  report. 
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In  PEP,  the  position  and  velocity  of  a  planet,  of  the  Earth -Moon  barycenter  or  of  the  Moon 
are  determined  as  functions  of  time  by  numerically  integrating  the  Encke  differential  equations 
for  these  quantities,  with  the  positions  of  perturbing  planets  being  determined  during  the  integra¬ 
tion  from  an  input  magnetic  tape.  The  partial  derivatives  of  position  and  velocity. with  respect  to 
masses  and  other  parameters  (not  initial  conditions)  are  determined  as  functions  of  time  by  nu¬ 
merically  integrating  the  differential  equations  for  these  quantities,  while  those  with  respect  to 
initial  conditions  are  determined  as  functions  of  time  either  by  assuming  that  they  are  equal  to 
the  partial  derivatives  of  position  and  velocity  with  respect  to  initial  conditions  in  the  elliptic  or¬ 
bit  osculating  to  the  true  orbit  of  the  body  at  the  initial  time,  or  by  numerically  integrating  the 
Encke  equations  for  these  quantities. 

We  feel  that  numerical  integration  on  an  electronic  computer  using  Encke 's  method  can  yield 
centuries  of  planetary  and  Earth-Moon  barycenter  ephemerides  of  the  accuracy  required  by  the 
observations  to  which  the  ephemerides  must  be  fitted.  In  the  case  of  the  Moon,  however,  nu¬ 
merical  integration  of  Encke' s  equations  might  not  yield  an  ephemens  accurate  for  centuries, 
although  it  certainly  v/ould  have  the  required  accuracy  for  decades.  Thus,  we  would  have  to 
manipulate  the  equations  further  into  a  form  that  would  yield  accurate  results  for  centuries  of 
numerical  integration! 

We  have  not  used  the  traditional  method  of  obtaining  planetary  motions  by  expansions  in  series 
for  a  number  of  reasons.  First,  the  higher  accuracies  we  require  necessitate  a  higher-order 
perturbation  theory  and  many  more  terms  in  the  truncated  series  than  v/ere  required  when  the 
present  ephemerides  were  generated.  Further,  the  equations  we  numerically  integrate  are  di¬ 
rectly  derivable  from  the  theory  of  gravitation  and  motion  employed,  while  many  operations  have 
to  be  carried  out  to  derive  the  series  solutions,  thus  introducing  the  possibility  of  error.  Finally, 
it  is  easy  to  introduce  additional  forces  in  the  numerically  integrated  equations,  whereas  the  con¬ 
sideration  of  the  effect  of  additional  forces  on  the  series  solutions  requires  much  effort. 

PEP  could  be  expanded  to  numerically  integrate  the  equations  for  the  motion  of  the  Earth 
about  its  center  of  mass,  in  addition  to  the  equations  for  the  motions  of  planets  around  the  Sun, 
of  . the  Earth-Moon  barycenter  around  the  Sun  and  of  the  Moon  around  the  Earth.  In  this  way,  all 
the  present  ephemerides  of  the  motions  in  the  solar  system  could  be  completely  and  rigorously 
revised,  instead  of  only  revising  the  ephemerides  of  the  motions  of  the  center  of  masses  of  the 
various  bodies,  assuming  the  present  expressions  for  the  rotation  and  precession-nutation  of  the 
Earth.  However,  even  if  we  assume  these  expressions,  significant  improvements  in  the  ephe¬ 
merides  of  the  motions  of  the  center  of  masses  can  be  made. 

Recent  radar  observations  at  Lincoln  Laboratory  and  elsewhere  of  Mercury,  Venus,  Mars 
and  the  Moon^  have  much  greater  accuracy  than  optical  observations  of  these  bodies.  However, 
optical  observations  have  the  advantage  of  having  been  made  over  a  period  of  several  centuries. 
Using  both  radar  and  optical  observations  to  improve  ephemerides  takes  account  of  the  stated 
advantages  of  both  kinds  of  observations.  In  addition,  the  dimensionality  of  the  space  of  observa¬ 
tions  is  increased  to  four  over  the  two  obtainable  using  only  optical  observations;  that  is,  radar 
observations  of  time  delay  and  doppler  shift  give  range  and  range-rate  measurements  in  addition 
to  the  two  angular  measurements  given  by  optical  observations. 

f  These  manipulations  have  since  been  performed  and  will  be  presented  in:  M.  E.  Ash,  "Generation  of  the  Lunar 
Ephemerii  on  an  Electronic  Computer,"  Technical  Report  400,  Lincoln  Laboratory,  M.l.T.  (24  August  1965). 

|The  Sun  and  perhaps  Jupiter  have  also  been  detected  with  radar,  but  the  results  of  such  observations  are  not  yet 
of  the  accuracy  or  nature  needed  in  improving  ephemerides. 
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With  PEP,  the  fact  that  fitting  ephemerides  to  observational  data  is  done  by  an  electronic 
computer  and  is  completely  automated  implies  that  more  accurate  ephemerides  will  be  generated 
than  if  traditional  methods  {largely  dependent  upon  hand  computation)  were  used,  even  with  ex¬ 
actly  the  same  observational  data  as  input.  Of  course,  more  observational  data  are  available 
now  than  when  the  present  ephemerides  were  generated,  even  without  counting  radar  observations. 

In  the  process  of  improving  ephemerides,  we  obtain  improved  values  of  the  various  param¬ 
eters  appearing  in  the  theory  {such  as  planetary  and  lunar  masses);  we  also  test  the  validity  of 
the  theory  employed.  Some  hypotheses  which  we  are  interested  in  testing  are: 

(1)  Does  the  Sun  have  a  detectable  second  harmonic  in  its  gravitational 
potential? 

(2)  Are  the  values  of  the  gravitational  constant  and  the  velocity  of  light 
functions  of  time  ? 

(3)  Is  there  an  advance  of  the  perihelion  of  the  orbit  of  Mercury  and  the 
other  planets  as  predicted  by  general  relativity? 

(4)  Is  fhe  general  relativistic  expression  for  the  time  delay  of  a  radar 
signal  passing  near  the  Sun  correct  ?* 

(5)  Can  atomic  time  be  identified  with  the  proper  time  of  general  relativity? 

The  prediction  of  the  advance  of  the  perihelion  of  Mercury's  orbit  has  been  verified  previously. 

Since  experimental  results  are  supposed  to  be  reproducible,  and  since  we  intend  to  check  this 

advance  with  more  accurate  data,  the  effort  we  make  to  do  this  is  not  without  merit.  In  order 

to  test  the  general  relativistic  effect  on  tlie  time  delay  of  the  radar  signal,  we  need  observations 

of  Venus  or  Mercury  at  superior  conjunction  at  the  frequency  of  Lincoln's  Haystack  radar 
9 

(8  x  10  cps).  These  have  not  yet  been  made,  but  hopefully  will  be  made  in  the  next  few  years. 


PRECEDING  PAGE  BLANK 


II.  ELLIPTIC  MOTION 


CHANGE  OF  COORDINATES 


In  P'ig.  1,  R  is  the  longitude  of  the  ascending  node  of  an  elliptic  orbit,  i  is  the  angle  of  in¬ 
clination  of  the  orbital  plane,  and  co  is  the  argument  of  perigee.  We  wish  to  find  the  transfor¬ 
mation  from  the  (x,  y,  z)  coordinate  system  to  the  (x,  y,  z)  coordinate  system,  whose  x-axis  is 
pointed  in  the  direction  of  perigee,  whose  y-axis  is  pointed  in  the  direction  of  motion  at  perigee, 
and  whose  z-axis  is  perpendicular  to  the  orbital  plane. 


Fig.  1.  Euler  angles. 


If  we  rotate  the  xy-plane  about  the  z-axis  through  the  angle  R,  we  obtain  an  (x',y\z’)  co¬ 
ordinate  system  related  to  the  original  one  by  the  equations 

x'  =  x  cos  R  +  y  sin  R 
y'  =  —x  sin  R  +  y  cos  R 
z'  =  z 

Rotating  about  the  x'-axis  through  the  angle  i,  we  obtain 


y"  =  y'  cos  i  +  z'  sin  i 
z"  =  — y'  sini  +  z'  cos  i 

Finally,  rotating  about  the  z"-axis  through  the  angle  w,  we  see  that 

x  =  x"  cos  w  +  y"  sinta 
y  =  — x"  sinui  +  y"  cosu 
z  =  z" 

The  net  result  of  these  three  transformations  is  seen  to  be 

x  =  (cos  R  cos  w  -  sin  R  sin  w  cos  i)  x  +  (sin  R  cos  w  +  cos  R  sin  w  cos  i)  y 
+  (sin  oj  sin  i)  z 

y  =  -(cos  R  sin  oj  +  sin  R  cos  w  cos  i)  x  +  (-sin  R  sin  w  +  cos  R  cos  cos  i)  y 
+  (cos  oj  sin  i)  z 

z  =  (sin  R  sin  i)  x  -  (cos  R  sin  i)  y  +  (cos  i)  z  .  (II-l) 

Equation  (II-l)  is  an  orthogonal  transformation,  so  its  inverse  transformation  is  given  by 
the  matrix  which  is  the  transpose  of  the  above  matrix.  Thus,  we  have 
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(II— 2) 


x  =  (cos  ft  cos  os  —  sin  ft  sin  03  cos  i)  x  —  (cos  ft  sin  03  +  sin  ft  cos  03  cos  i)  y 
+  (sin  ft  sin  i)  z 

y  =  (sin  ft  cos  03  +  cos  ft  sin  03  cos  i)  x  +  (—sin  ft  sin,  03  +  cos  ft  cos  03  cos  i)  y 
-  (cos  ft  sini)  z 

z  =  (sin  03  sin  i)  x  +  (cos  03  sin  i)  y  +  (cos  i)  z 
These  formulas  agree  with  those  given  in  Ref.  2. 


B.  DETERMINATION  OF  POSITION  AND  VELOCITY  AT  A  GIVEN  TIME 
FROM  ORBITAL  ELEMENTS 


The  differential  equation  system 


d2y2  =  _  Hyj 
dt2  T3" 


j  =  1,2.3 


(II— 3) 


/  j  2  2  2  3  2 

where  (i  >  0  and  p  =  v  (yap  +  (y  )  +  (yj  ,  represents  the  motion  of  a  body  traveling  in  a  conic 
section.  We  shall  assume  that  this  conic  section  is  an  ellipse. 

The  various  parameters  used  to  describe  elliptic  motion  are: 


a  =  semimajor  axis 
e  =  eccentricity  (0  ^  e  <  1) 
i  =  orbital  inclination  (0^  i<  180°) 
ft  =  longitude  of  the  ascending  node  (0  ^  ft  <  360°) 

03  =  argument  of  perigee,  measured  along  the  orbital  plane  from  the  ascending 
node  (0°  ^  03  <  360°).  (Note:  the  longitude  of  perigee  is  the  quantity 

OJ  =  SI  +  03.) 


iQ  =  initial  mean  anomaly  at  time  tQ 
£  =  mean  anomaly  at  time  t 
f  =  true  anomaly  at  time  t 
u  =  eccentric  anomaly  at  time  t 
n  =  mean  motion 
p  =  semilatus  rectum 


(0°  ^  t  <  360°) 
(0°  ^  £  <  360°) 
(0°  <  f  <  360°) 
(0°  u  <  360°) 


The  quantities  (a,  e,  i,  ft,  03,  £  )  are  the  elements  of  the  elliptic  orbit  and,  along  with  the  time, 

0  1  2  3  •  1  •  2  . 3 

completely  determine  the  position  and  velocity  (y  ,  y  ,  y  ,  y  ,  y  ,  y  )  in  the  orbit. 

Let  (x,  y)  be  the  Cartesian  coordinate  system  in  the  orbital  plane  whose  x-axis  points  in  the 
direction  of  perigee  and  whose  y-axis  points  in  the  direction  of  motion  at  perigee.  The  following 
equations  are  derived  in  Ref.  3  (the  formula  numbers  in  square  brackets  are  those  in  the  reference). 


(y1)2  +  <y2)2  +  (y3)2  =  n(|-  |)  [26,41] 


(II— 4) 


=  a(l  -  e2) 

(II  -  3) 

P 

1  +  e  cosf 

[44] 

(II-6) 

=  a(l  -  e  cos u) 

[51] 

(II—  7) 

=  a'3/2 

[43] 

(II— 8) 

=  'o  +  n|t-V 

[50] 

(II— 9) 
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l  -  u  -  e  sin  u 

(49] 

(II— 10) 

f  / 1  +  e  ,  u 

tari-r  =  - -  tan -7 

2  v  1  -  e  2 

[52] 

(II-ll) 

x  =  a(cos  u  —  e) 

[54] 

(11-12) 

/,  2  . 

y  =  a  v  1  -  e  sin  u 

[55]  . 

(11-13) 

These  are  the  basic  formulas  which  we 

(II- 10)  with  respect  to  time,  we  find  that 

shall  assume  to  be  known. 

Differentiating  (II-9)  and 

d£  .  du 

n=  dt  =  (1-ec°su)  dT 

Thus, 

du  _  na 
dt  ”  p 

(11-14) 

Equations  (11-12),  (11-13)  and  (11-14)  together  give 

2 

—  na 

x  - - sinu 

P 

(11-15) 

2  /,  2 
—  na  v  1  —  e 
y  =  -  cos  u 

n 

(11-16) 

1  2  3-1  2-3 

Suppose  we  are  given  (a,  e,  i,  ft,  w,  £q),  and  we  wish  to  find  (y  ,  y  ,y  ,y  ,y  ,y  )  at  time  t. 
Equation  (II-9)  determines  the  mean  anomaly  £  at  time  t,  from  which  Kepler's  equation  (II - 10) 
(solved  by  iteration)  gives  the  eccentric  anomaly  u  at  time  t.  Then,  formulas  (II- 7),  (11-12), 
(11-13),  (11-15)  and  (II- 16)  determine  x,  y,  ij,  j  at  time  t.  If  we  define 

=  cos  ft  cos  w  -  sin  ft  sin  w  cos  i 

b12  =  —cos  ft  sinw  -  sin  ft  cosw  cosi 

b0.  =  sin  ft  cos  w  +  cos  ft  sinw  cosi 

b22  =  —sin  ft  sin  w  +  cos  ft  cos  w  cos  i 

bOJ  =  sinw  sini 
31 

b32  =  cosw  sini  (11-17) 


we  finally  obtain,  by  (II-2), 


^  ’  V  * +  bj2  v 

y>‘bni*b^. 


j  =  1,2,3 


(11-18) 


C.  DETERMINATION  OF  ORBITAL  ELEMENTS  FROM  POSITION 
AND  VELOCITY  AT  A  GIVEN  TIME 

1  2  3  •  1  •  2  -3 

Suppose  we  are  given  (y  ,  y  ,  y  ,  y  ,  y  ,  y  )  at  time  t,  and  are  required  to  find  (a,  e,  i,  R,  w,  £  ), 
with  £q  being  the  mean  anomaly  at  time  t  . 


i 

i  1 
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We  define 


2  ,  1.2  J  .  2.2  ,  .  3,2 
p  -  (y  )  +  (y  )  +  (y  ) 

2  .-1,2  ,  ..2,2  ,  ..3,2 

V  =  <y  )  +  (y  )  +  (y  ) 

-  -  1.1  ,  2.2  ,  3.3 

P  •  V  =  y  y  +yy  +  y  y 

Then  the  vis  viva  integral  (II-4)  immediately  determines  one  of  the  elements: 


(11-19) 


a  "  2,i  2 

— - v 

P 

Differentiating  (II-7)  with  respect  to  time,  we  obtain 


(11-20) 


dp  .  du  na  e  smu 

-sf  ae  sinu  ~rr  =  - 

dt  dt  p 


by  (11-14).  Since 


dp  _  ~p  •  v 
3t  -  ’  p 


this  gives 


P  .  V 

e  sxnu  =  - — 


(11-21) 


Further,  (II-7)  can  be  put  in  the  form 


e  cos  u  =  1  -  i- 
a 


(11-22) 


The  simultaneous  solution  of  (11-21)  and  (11-22)  determines  e  and  u.  The  mean  anomaly  at 
time  tQ  is  then  found  from  Kepler's  equation  (II- 10),  and  formulas  (11-12)  through  (II- 16)  deter¬ 
mine  x,y,x,y.  Solving  (11-18)  for  the  b^k,  we  have 


(xy  -  xy) 


(y3y  -yJy) 


(xy  -  xy) 


(y 3 x  -  y3x) 


j=  1,2,3 


which  can  be  put  in  the  form 


bji  =  yj(c-fH,-yj(^) 

,  i  /  sinu  \,  -j  /  cos  u  -  e  \ 

j2  ' 


j  =  1,  2,  3  . 


(11-23) 


These  formulas  agree  with  those  in  Ref.  4,  in  which  \TjI  t  was  used  as  the  time  variable  instead  of  t. 

By  means  of  the  above  methods  we  can  thus  determine  a,  e,  £0,  and  b^t.  (j  =  1,  2,  3;  k  =  1,  2) 
given  the  position  and  velocity  (y  ,y  ,y  ,y  ,  y  ,y  )at  time  t.  From  these  quantities,  we  can  then 

1  2  3-1  2.3 

determine  the  values  of(y  ,y  ,y  ,y  ,y  ,y  )  at  any  time  from  the  formulas  given  inSec.  II-B  above. 
To  determine  i,  B  and  w,  we  must  solve  (11-17).  First,  we  have 


8 


0°  i-$  180 


(11-24) 


sini  =  J 


,2  x  K2 

b31  b32 


COSi  =  bllb22-bi2b2lJ 


Further,  if  i  ^  0°  or  180°,  we  see  that  co  is  determined  by  the  relations 


sin  co 


cos  co 


31 


J 


v  2  ,  .2 

b31  +  b32 


}  0°  4  w  <  360c 


r* 


uzz 

~ 7~T 

31  +  b32 


(11-25) 


Finally,  we  see  that  ft  is  determined  by  the  relations 


sin  SI 

=  ^>2i  cos  w  — 

b^2  sin  co 

0°  ^  ft  <  360° 

cos  SI 

=  b,  .  cos  co  — 
11 

b^2  sin  co  j 

180% 

the  equations 

in  (11-17)  assume  the  form 

bll 

=  cos  (ft  ±  co) 

b^2  =  ±  cos  (ft  ±  co) 

b12 

=  —sin  (co  ±  ft) 

b3l=0 

b21 

=  sin  (ft  ±  co) 

b32  =  0 

(11-26) 


where  the  plus  sign  of  the  ±  symbol  is  to  be  used  if  i  =  0°,  and  the  minus  sign  is  to  be  used  if 
i  =  180°.  Thus,  ft  and  co  are  indeterminate  when  i  =  0°  or  180°.  We  might  make  the  conven¬ 
tion  that  when  i  =  0°  or  180°,  we  set  ft  =  0  and  determine  co  by  the  relations 

cos  co  =  b , 


'll 
sin  co  =  -b 


12 


(11-27) 


D.  PARTIAL  DERIVATIVES  OF  POSITION  AND  VELOCITY 
WITH  RESPECT  TO  ORBITAL  ELEMENTS 

12  3  12  3 

Regarding  the  position  and  velocity  (y  ,  y  ,  y  ,  y  ,  y  ,  y  )  at  time  t  as  functions  of  the  orbital 
elements  (a,  e,  i,  ft,  co,  (Q),  v.e  derive  the  following  equations. 


9yJ  _yj_3y]  ,  . 

da  ~a  2  a  o' 


ar 

8a 


1  3ex(t  -  tQ) 
2a  +  2ap 


♦  L,J 

2 P  \  Jl 


■TTez 


—  b.7J  1  -  e2  x 

l2 


j  =  l,  2,  3  (11-28) 


^  =  — ab. .  - 

9e  31  l  - e 


_  b  v+  y3  sinu 

2  j2  y  +  n 


.  .  i  .i  /  x.  .  x  b..ax  cosu  /x._  \ 

iT  =  Z1  fa  cosu  4  exjmu)  +  jj -  +  b  / H  -  ) 

9e  p  \  n  )  p  J2  \  P  i  _  e2/ 


j  =  1,  2,  3 


(11-29) 


..  '  : 


1 


3y3  _  yjJ. 

9fo  "  n 


r  3  \r  J  Ck  \r 


3yJ  _  yJex  A  a  /  u 
at-  -  —  +  p(-bji 


+  b.,  v  1  -  e2  x 
J2 


j  -  1,2,3 


(11-30) 


%U=(sin£2)  y3 


=  (sin  £2)  y 


=  -(cos  £2)  y3 


=  -(cos  £2)  y3 


(H-31) 


=  (sin  w  cos  i)  x 
3i 

+  (cos  w  cos  i)  y 


LC  =  (sinw  cos  i)  x 
3i 

+  (cos  w  cos  i)  y 


iiL  =_y2 

3  £2  y 


3y  1 

3ir  =  y 


*  • 

9y  =  ,-,i 


(II-32) 


iy3  -  o 

ll£2  _  0 


ay3  _  o 

W  "  0 


9yJ  _  ,  _ 

T&  ‘VV 
2£i  -b.,s-b..y 

do)  ]2  jlJ 


j  =  1.  2.  3  . 


(11-33) 


To  show  the  validity  of  the  above  formulas,  we  note  that  in  (11-18)  the  b^k  are  functions  of 
i,  £2,  w,  and  x,  y,  x,  y  are  functions  of  a,  e,  lQ,  t.  Therefore,  we  first  differentiate  (11-17)  with 
respect  to  i,  £2,  w,  and  find  that 

abll 

—  =  sin  £2  sinw  sini  =  (sin  £2)  b,  , 

3i  '31 


sin  £2  cosw  sini  =  (sin  £2)  b. 


-gj-  =  —cos  £2  sin  w  sini  =  -(cos  £2)  b^ 


—  =  -cos  £2  cos  w  sini  =  —(cos  £2)  b,9 
91  oc 


=  sin  w  cos  i 


=  COS  O)  cos  1 


(11-34) 


-sin  ft  cos  co  —  cos  ft  sin  co  cos  i  =  — b. 


g  ^  =  sin  ft  sin  co  —  cos  ft  cos  co  cos  i  =  -b ^ 


Tn"  =  cos  ft  cos  co  —  sin  ft  sin  co  cos  i  =  b .  . 
d  ft  11 


~dQ~  =  —cos  ft  sin  co  —  sin  ft  cos  co  cos  i  =  b^ 


=  0  ; 


(11-35) 


— —  =  —cos  ft  sin  co  -  sin  ft  cos  co  cos  i  =  b,0 
8co  12 


—  =  —cos  ft  cos  co  +  sin  ft  sin  co  cos  i  =  — b .  . 

9co  11 


— —  =  —sin  ft  sin  co  +  cos  ft  cos  co  cos  i  =  b~- 
3co  22 


-r —  =  -sin  ft  cos  co  -  cos  ft  sinco  cosi  =  —  b,, 
9  co  21 


— r —  =  cos  co  sini  =  b,-> 
9  co  32 


— —  =  —sinco  sini  =  -b,  . 
3<o  31 


(11-36) 


Then,  (11-34)  through  (11-36)  and  (11-18)  together  show  that  formulas  (11-31)  through  (11-33)  are 
valid. 

Next,  by  (II-8),  we  have 


9n  _  _3  n 
3a  2  a 


(11-37) 


Further,  (II-9)  and  (II- 10)  give 


_2  n  ,  v  =  9«  =  £u  9u 

2  a  U  o  3a  3a  3a 


0  =  jp  =  —  sin  u  -  e  cos  u 

3e  3e  9e 


,  3  (  3 u  3u 

1=3£  ~  dt  e  cos  u  g  i 
o  o  o 


By  (II-7),  these  equations  imply  that 


11 


JL 


-I 


ft-. 


■  2 


9u 

9a 


3  n(t  “  to) 


2  p 

9u  _  a  sinu 
9e  p 


9u  a_ 

3fo"  P 


(11-38) 


From  (II— 7)  and  (11-38),  it  follo.ws  that 

.  .  -  nae(t  - t  )  sinu 

9 p  ,,  .  ,  .9up3  o 

— t-  =  (1  —  e  cos  u)  +  ae  sin  u  —  =  —  —  or - 

9a  9a  a  2  p 


9  p  _  ,  9u  ,  a  e  sin  u 

—  =  —a  cosu  +  ae  sinu  o—  =  — a  cosu  +  - 


9e 
3  P  . 

9t: 


ae  sinu 


9e 


9u  a  e  sin  u 


9 1 


(11-39) 


Having  obtained  the  above  formulas,  we  use  (II- 12)  to  (II- 16)  to  find  the  derivatives  of  x,  y, 
x,  y  with  respect  to  a,  e,  l  . 

9x  ,  .  9u  x  3  x  . 

9a  9a  a  2  a  o 

9x  _  _  9u  _  _  ,  xsinu 

t—  =  —a  —  a  sin  u  —  =  —a  +  - 

9e  9e  n 


9x  _  .  9u  x 

9 1  a  sinu  9£  n 

o  o 


(11-40) 


=  Jl  -  e2  sin  u  +  a  J 1  -  e2  cos  uljs|-l|(t-t  ) 
33  3  3  3  Z  3  O 


9y  ae  sinu  .  j. 2 
— -  - -  +  a  N  1  -  e  cos  u 


9e 


9u 

9e 


ey 


1  -  e 


y  sinu 
n 


9y  f.  2  9u  y 

rf  =  aNll-e  cosu  TT7- =  4  : 


9£ 


9£  n 
o 


(11-41) 


•  2  2  2 
3x  _  _  a  sin  u  9n  _  2na  sinu  _  na  cos  u  9u  na  sinu  9 p 

9a  ”  p  9a  p  p  9a  2  9a 

P 

2  2 

,  ,  3n  a  cos  u  (t  -  t  )  ,  r .  3nae  sin  u  (t  —  t  ) 

—  ,3.2,.  o—ll  o 

=  X  (—  5 —  +  — )  +  - -  —X  I - - ~ - 

2a  a  0.2  a 


2P 


2p 


1  ,  3ex(t  ~  l0) 


2a  + 


2a  p 


3ny(t  -  tQ) 


2p  v  1  -  e 


(11-42) 
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2  2 
9x  _  -na  cosu  9u  na  sin 

9e  p  9e  2 

P 


u  dp  _  ax  cosu  ,  x  /  _  exsinu\ 

9e  p  p  \  n  / 


—  2  2  ^  ju 

9x  _  na  cos u  8u  A  na  sinu  3 p  ay  x  xex 

3,o"  0  8,o  p2  no’  T. - 5  np 

^  p  v  1  -  e 


(11-43) 


9y  _  a 
9a 


“  *fi  —  e2  cos  u  9n  2na  ’J i  —  cos  u  _  na^  1  —  sir,  u  9u_ 


2  / 2 

_  na  v  1  —  e  cos  u  9p 

2  9a 

P 


2  2/  2 

,  ,  3n  a  vl-e  sinu(t  —  t  )  ,  r.  3naesinu(t  — t  ) 

=  n-At||  +  - - - 2.-7  |± - 5 - 2- 


2p2  '  1* 


.  .  3ex(t  —  t  )]  3n  n  1  —  e  x(t  -  t  ) 

=  —v  —  +  - —  —  _ 2_  : 

y  1 2a  2a p  2 p 


(11-44) 


2  2  r  2  .  2  /  2 

i  e  cos  u  na  n  1  —  e  sm  u9u  na  v  1  —  e 


9y  _ _  na  e  cosu  _  na  v  1  —  e  sinu  9u  _  na  n  1  —  e  cos u  9 p 

3e  / - T  o  9e  2  'We 

pj  l-eZ  ‘  P 


=  --JLy...  +  211  +  y  (a  COSU+  £2Lsmu} 
U-e2)  0  0 


—  2  /  2  2  / 2  —  —  / 2  — 
9y  na  V 1  -  e  cosu  9p  na  v  1  —  e  sinu  9u  _  yex  ,  a  V  1  —  e  x 

9  £  "  2  9£  p  8«  "  np  p 

O  p  o  r  o  r  r 


(11-45) 


Equations  (11-40)  through  (11-45)  and  (11-18)  together  show  that  formulas  (11-28)  through  (11-30) 
are  valid. 

E.  PARTIAL  DERIVATIVES  OF  ORBITAL  ELEMENTS  WITH  RESPECT 
TO  POSITION  AND  VELOCITY 

Regarding  the  orbital  elements  (a,  e,  i,  SI,  w,  £  )  at  time  t  as  functions  of  the  position  and 

123.123  °  0 

velocity  (y  ,y  ,y  ,  y  ,  y  ,  y  )  at  time  t,  we  derive  the  following  equations. 


9a 

->  2  i 
_  2a  yJ 

ay3 

3 

P 

9a 

,  2.  i 
_  2a  yJ 

ay3 

P 

9e 

sinu 

9yj 

na 

9e 

_  sinu 

ay3 

na 

j  =  1.  2,  3 


(11-46) 


9e  _  sinu  y3  _  (p  ■  v)  y3  +  y3  cosu  .2  _  _1. 

3y  j  na  a  3  J  p  (P  a' 


j  =  1,2,3 


(11-47) 
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I 


!!o 

9y3 


cos  u  —  e 
nae 


y3  (p  •  v)  yJl  y3  sinu  ,2  1,  , 

"a  — p  (p~a)  + 


3nay3(t  -  tQ) 


3£ 

ay 


o  _  cos  u  -  e 

j 


nae 


y3  _  (p  •  v)y3l  _  2py 3  sinu  s 
a  (i 


3nay  3(t  -  tQ) 


3=1 


— •  =  ( 
ay3  ' 


9b 


31 


3b 


9y 

b 


3 


sin  a;  + 


32 


9yJ 


COS  CO 


j  cos  i  +  ^ 


h  9b12 

b12  _j  +  b21 


9yJ 


3yJ 


9b,,  3b, 

_  u  _ ; 

11  n  3  22  . 

9yJ  3y 


yj  sini  ,  j  =  1, ....  6 


where  (y1,  .  . .  ,  yb)  =  (y  , . . .  ,y  ),  and  where  the  3b^/3y3  are  given  in  (11-58)  and  (II 
The  same  remarks  apply  to  the  following  two  sets  of  equations. 

9co  1 


9y3  smi  \9y3 


9b31  9b32  . 

coso) - —  sinco 


ay 


3 


3  0 
9y3 


f3b 


21 


3b 


9y ' 


COS  0)  - 


22 


j  =  1, .  .  . ,  6 
3b, 


3y 3 


sin  w  )  cos  0  — 


'11 


,3 


COS  0) 


9b 


9y 


12  .  \  .  _  3  co 

T  smcoi  sinfi  —  r—  cosi 
3  /  9a 


3yJ 
j  =  1,  ....  6 


Equations  fII-46)  follow  directly  from  formula  (11-20).  To  derive  (11-47)  and  (II- 
first  differentiate  (11-21)  and  (11-22),  obtaining 


9(e  sinu)  _  y3  (p  •  v)  y3 


9y- 


na 


nap 


3(e  sinu)  _  y3  (p  •  v)  y3 
-  1  ■  *“  — ^  ““  ■  —  — - 


3  =  1>  2,  3 


3y 


3 


na 


nap 


3(e  cosu)  2y3  y3 

-7- 

3(e  cosu)  _  2py3 


9y 


3 


3  =  1.2,3 


For  any  variable  a,  we  have 

9e 


3(e  sinu) 
da 


9u 


=  sinu  ~  +  a  cos  u  tt— 
9a  3a 


9(e-cPiu)  =  cosu  I®  -e  sinu 
3a  3a  3a 


which  can  be  put  in  the  form 

=  sinu  9(esinu)  9(ecosu) 

3a  3a  3a 


U  =  J  fcosu  -  sinu 


9(e  sinu) 
9a 


9(e  cosu), 
9a 


,2,3  (11-48) 

(11-49) 

59)  below. 

(11-50) 

(II-51) 
-48),  we 

(11-52) 

(H-53) 


(11-54) 
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Equations  (11-52)  through  (11-54)  together  imply  the  validity  of  formulas  (11-47)  and  of  the  following 
formulas. 


3  u 

9y3 


cos  1 
nae 


u  fyJ  _  (p  •  v)  y3 


9u  _  cos  u 

„  •  j  nae 
3yJ 


y3  _  (p  •  v)  y3 
a  P 


y3  sinu  ,2  _  1. 
ep  'p  a; 


_  2py3  sinu 

ep 


j  =  1,2,3 


(H-55) 


Let  a  be  any  variable.  Differentiating  (II— 9)  a.,d  (II- 10)  with  respect  to  a,  we  see  that 


3  £ 

_ o  3n_  ,,  _  .  ,  d±  _  3u.  _  3(e  sinu) 

3a  ace  o-aa_3o;  3a 


(11-56) 


where,  by  (II-8), 


3n 

da 


3  n  3a 
2  a  da 


(H-57) 


Formulas  (11-48)  then  follow  from  (11-52)  and  (11-55)  to  (11-57). 

Next,  we  differentiate  (11-23)  with  respect  to  y3  and  y3,  obtaining 

3b  k  i  • k 

kl  _  cosu  ^  _  y  y J  cos  u  _  y  sin u  3a 

3y3  P  kj  p  2na  ay3 


k  •  k. 

/  y  sin  u  x  y  cos 

\  P 
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u\  3u 

i  3yj 


9bkl  _  _  sinu  ^  _  yk  sinu  3a  _  / yk  sinu  A  yk  cosu\  9u 

ay3  "  na  k3  2na2  3y3  V  P  na  >  dy\ 


ab 
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'sinu 


k  j  . 

J  V  J  <21 


ay 


•fi  -  e2 


kj 


y  yJ  sinu  ,  y  (cos u  -  e)  3a 

—  - - -  T  - T -  - - 


2na ' 


3yJ 


(k  k  \ 

y  cosu  _  y  sinul  3u 

p  na  /  ay3 


yke  sinu  ,  y^  [e(cos u  -  e)  _  . 

2.  na  ..  2. 

p(l-e)  l(l-e) 


3b 


k2 


(1  -  e  ) 


•  k. 


3e 

3y3 


j  =  1,2,3 
k  =  1,  2,  3 


(11-58) 


3y3 


1 _  /(cos  u  -  e)  g  +  y  (cos  u  -  e)  d  a 


>  j  =  1,2,3 
f  k=  1,2,3 


Ji  -e2 


na 


kj 


2na 


3yJ 


k  •  k 

+  |  j  cos  u  _  y  sin 
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na 


u\  au 
'  3y3 


yke  sinu  ,  y^  [e(cos u  -  e)  _  , 

“77  “Z7  na  2~ 

p(l  -  e  )  l  (1  -  e  ) 


Here,  the  Kronecker  delta  <5^  is  defined  by 

[0  ifj^k 


3e 

3y3 , 


(11-59) 


kj 


1  if  j  =  k 


(11-60) 
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Differentiating(II-24)  with  respect  to  y3  (j  =  1, . . . ,  6),  we  see  that 


cos i 


3y3  dy3 


sinw  + 


.  .  di  u  9b22  ,  ,  9bU  ,  9b21  ,  9b12 

-smi - -  =  b.  .  - —  +  b00  - . - b ^ - b, ,  - t- 

ay3  11  ay3  22  ay3  12  ay3  1  ay3 


(11-61) 


Multiplying  the  first  equation  of  (II-61)  by  cos  i  and  the  second  equation  by  sini,  and  subtracting, 
we  obtain  (11-49).  Further,  differentiating  either  the  first  or  the  second  equation  of  (11-25) 
yields  (11-50).  Finally,  by  differentiating  (11-26)  we  see  that 

an  abo-i  ab?2  0m 

cos  ft  - r  = - 1- -  cos  u - s*na)  —  (^21  s*na)  +  ^22  cos  - r 

3y 3  9y3  3y3  3y3 


««  3b  i  i  3b  a  a. 

-sin ft  -2— =•  =  - :=■  cosw - ^  sinw  -  (b, ,  sina>  +  b.?  cosw)  — r 

ay3  ay3  ay3  11  w  ay3 


(11-62) 


Multiplying  the  first  of  these  equations  by  cos  ft  and  the  second  by  sin  ft,  and  subtracting,  we 
obtain  (11-51).  Here,  we  use  the  fact  that  by  (11-17) 

cos  i  =  (b21  sinw  +  b22  cos  u)  cos  ft  —  (b^  sinu  +  b^z  cos  w)  sin  ft 

F.  CHECK  OF  ELLIPTIC  ORBIT  FORMULAS 

The  elliptic  orbit  formulas  derived  in  the  preceding  sections  were  used  to  write  computer 
subroutines  needed  in  Encke's  method  for  the  numerical  integration  of  the  equations  of  motion 
and  the  equations  for  the  partial  derivatives  with  respect  to  initial  conditions  of  a  planet.  These 
computer  subroutines  then  enabled  us  to  check  the  validity  of  the  elliptic  orbit  formulas  in  the 
following  manner. 

First,  note  that  the  position  and  velocity  (y1,  .  .  . ,  yS  in  an  elliptic  orbit  satisfy  the  dif¬ 
ferential  equations  system 

dyk  k+3 

-ar  =  y 


.  k+3  k 

oy  _  _  ny 

~w  V 

k  k 

y  =y„ 


k  =  1,2,  3 


(11-63) 


k+3  k+3 


when  t  =  t 


Let  (0 \  . . .  ,0^)  denote  the  elements  (a,  e,  i,ft,  w,  t  )  of  the  elliptic  orbit.  Differentiating  system 
(11-63)  with  respect  to  /33,  we  obtain 

d(8yk/a p3)  _  3yk+3  ] 


dOyk+3/afl3)  _  p  h 

ar~  '  i  " 

p  \  i 


Z  y1  — j  -  *4 
P  £=1  w3  w} 


k  =  1,  2,  3 

j«l . 6  .  (11-64) 


a/33  a/33 


ayk+3 


when  t  =  t 
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L,-  .  . 


k  k  /  i 

Given  the  orbital  elements,  we  calculated  numerically  the  quantities  y  and  9y  /9/3J  (j,  k  =  1, . . . ,  6) 

l’ 

•  for  selected  values  of  time  using  the  formulas  of  Secs.II-B  and  II-D  above.  Also  using  these  for- 

,  mulas  to  determine  the  initial  conditions  yQ  and  dyo/dPJ  (j,  k  =  1, . . . ,  6),  we  numerically  inte¬ 

grated  the  42  differential  equations  (11-63)  and  (11-64).  The  results  of  the  numerical  integration 

*  and  the  calculations  from  the  formulas  agreed  to  within  the  accuracy  expected  of  the  numerical 

integration,  showing  the  validity  of  the  formulas  in  Secs.II-B  and  II-D.  This  was  also  a  check 

k  k  /  i 

of  the  computer  subroutines  to  calculate  y  and  9y  / 9/3 J  (j,  k  =  1, . . . ,  6),  and  of  the  numerical 
integration  subroutine.  (We  had  to  write  our  own  numerical  integration  subroutine,  since  there 
was  none  available  which  used  double  precision  arithmetic  operations.) 

Next,  given  the  orbital  elements,  we  calculated  the  position  and  velocity  at  the  initial  time 
from  the  formulas  in  Sec.  II-B.  Then,  using  the  formulas  of  Sec.  II-C,  we  calculated  the  orbital 
elements  from  this  position  and  velocity  and  observed  that  the  final  and  starting  orbital  elements 
agreed  to  witnin  the  roundoff  error  expected  in  the  calculations. 

Finally,  given  the  orbital  elements,  we  calculated  the  matrix  (dy^/dpi)  from  the  formulas 
in  Sec.  II-D  and  the  matrix  (9/^/9^)  from  the  formulas  in  Sec.II-E  for  selected  values  of  time. 
Multiplying  these  two  matrices,  we  found  that  the  result  differed  from  the  identity  matrix  to 
within  the  roundoff  error  expected  in  the  calculations. 

All  the  elliptic  orbit  formulas  derived  above  are  valid  for  all  choices  of  the  orbital  elements 
(a,  e,  i,  Q,  ),  except  that  son. ~-of  the  formulas  in  Sec.  II-E  for  the  matrix  (dp^/dy^)  become 
indeterminant  for  e  =  0  or  i  -  0°,  180°  because  here  the  matrix  (9y ^/dp^)  is  singular.  Now,  the 
only  time  that  the  formulas  for  (dp^/dy^)  a  re  used  is  when  we  are  integrating  the  equations  for 
the  partial  derivatives  of  the  position  and  velocity  of  a  planet  with  respect  to  initial  conditions 
using  Encke's  method,  and  we  wish  to  change  Encke  orbits.  It  is  not  very  probable  that,  at  the 
instant  when  the  Encke  orbit  is  changed,  the  orbital  elements  would  be  right  at  one  of  the  critical 
%  points.  In  the  case  of  the  integration  of  the  planetary  motions,  this  could  never  happen. 


f 
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III.  NEWTONIAN  EQUATIONS  OF  MOTION 

A.  TWO-BODY  PROBLEM  (PLANET-SUN)  WITH  PERTURBING  FORCES 

Uc  wish  to  find  the  equations  of  motion  of  a  celestial  body  about  a  central  body,  with  this 

motion  being  perturbed  by  N  other  bodies  {whose  positions  are  known  as  functions  of  time)  and 

by  forces  dependent  on  the  position  and  velocity  of  the  given  body  relative  to  the  central  body. 

Let  the  subscript  s  denote  the  central  body  (s  =  Sun),  p  the  given  body  (p  =  planet),  and 

j  (1  <  j  <  N)  the  Jth  perturbing  body.  Lot  y  denote  the  gravitational  constant,  and  suppose  that 
1  Z  3 

(x  .  x  ,  x  )  is  an  inertial  coordinate  system.  We  make  the  following  notational  conventions: 
xgk  =  coordinate  of  central  body,  etc. 

k  k  l  j. 

x.  =  x  —  x  =  coordinate  of  j  relative  to  s,  so  that  x  =  — x  etc. 

Jfa  J  s  Js  sj 

1  sj  ~  rjs  =  dis^ance  between  s  and  j,  etc. 

F  =  ktk  coordinate  of  additional  force  on  p 

P  1 

M  -  mass  of  s,  etc. 

Then,  by  Newton's  laws  of  motion  and  gravity,  we  have 


J-  --  y™sMP'-f  +VMS  Z 

rps  j  =  l  rjs 


k  =  1,  2,  3 


(III-l) 


U  ,X  X  X .  . 

M  — =  yM  M  -4E  +  yM  Y  M .  +  F 

p  dt2  -  p  7^  p .  ,  j  T7  p 

ps  J=1  JP 

By  dividing  the  above  equations  by  M  and  M  ,  respectively,  and  bv  subtracting  the  first  equa- 

s  p 

tion  from  the  second  equation,  we  find  that  the  equations  of  motion  of  p  relative  to  s  are 


d2xk 

- ?-  =  ~VM 


(-  *  ^  *  °k  *  w~  Fpk 


k  =  1,  2,  3 


(III  —  2) 


where  the  perturbing  planet  term  is  given  by 

N  /  k  k  \ 

,  M .  /  x .  x .  \ 

«k  =  ^s  2)  i^HF  --J? 

.  .  s  \  r.  r.  y 

J=1  \  JP  js/ 


k  =  1,  2,  3 


( I II  —  3 ) 


Kffccts  contained  in, the  (l/M  )  F*  term  in  (III-2)  determined  in  this  report  are 

Kk  =  general  relativity  effect  [see  (IV -52) } 

S  =  second  harmonic  of  the  Sun  [sec  (III-50)]  , 
u 

If  we  let  F^  denote  forces  acting  on  the  planet  in  addition  to  those  enumerated,  (1II-2)  becomes 


4-  =■  -VMS  (l  +  +  »k  +  Sk  +  ■£-  Fk  .  k  =  1.  2,  3 

'  s '  r  p  * 


(III— 1) 
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B.  THREE-BODY  PROBLEM  (EARTH-MOON-SUN)  WITH  PERTURBING  FORCES 

We  wish  to  find  the  equations  of  motion  of  the  Earth-Moon  barycenter  about  the  Sun  and  of 
the  Moon  about  the  Earth,  with  these  motions  being  perturbed  by  N  planets  (whose  positions  are 
known  as  functions  of  time)  and  by  additional  forces  Fg  and  Fm  acting  on  the  Earth  and  Moon, 
respectively.  We  assume  that  these  forces  are  expressible  in  terms  of  the  relative  positions 
and  velocities  of  the  Sun,  the  Earth  and  the  Moon. 

Let  the  subscript  s  denote  the  Sun,  e  the  Earth,  m  the  Moon,  and  j  (1  <  j  <  N)  the  jth  per¬ 
turbing  planet.  Otherwise,  the  notation  in  this  section  is  the  same  as  in  Sec.  III-A  above.  New¬ 
ton's  laws  of  motion  and  gravity  then  give 


.2  k  k 

d  x  x 

— =  yM  ~ -  +  yM 

dt2  e  r3  m 

es 


N 


J  y. 


ms 


j=*  "js 


.2  k 
d  x 
_ e_ 

dt2 


,2  k 
d  x 

m 


x 


=  yMg  -SS  +  yM 
res 


,k  N  k 

+  y  V  M.  J®.  +  -ji-  F1 
m3  '  Lj  i  3  M  e 

rme 


je 


dt 


=  yM 


s¥^Me 
r 

ms 


j=1 

k  N  k 

x  _  x. 

em 


r +y  l  MiJr-  +  w~ F 

J  r*  m 

j=1 


me 


r . 


m 


k  =  1,  2,  3 


(HI— 5) 


Let  the  subscript  c  refer  to  the  center  of  mass  of  the  Earth-Moon  system.  Thus, 


M. 


M  +  M 
e  m 


,  M  .  M 
k _ e  k  _ m 

xc  “  M  xe  +  M 
c  c 


m 


es 


xk  ^ 
os  M 

c 


me 


k  k  k 

X  K  =  x  K  +  — —  X  X 

ms  ‘  cs  M  me 
c 


k  =  1,  2,  3 


(III-6) 


By  multiplying  the  second  equation  of  (III-5)  by  M  /M  and  the  third  equation  by  M  /M  ,  and  by 

c  c  m  ^ 


adding  the  results,  we  obtain 

M 


,2  k 
d  xc 


M 


N 


dt 


=  - S  +  yM  - 21  .  Sm  +  y  V  M  _ 6 

yi  s  M„  3  yi  s  M  3  y  Li  lvli  M 


c  r 


es 


c  r 


ms 


Mm  xim  \  1  (Fk  „k 

+  T5T  r  )  +  m~(F«  +  F™> 


m_  x.1: 


c  r. 
jm/ 


m' 


(HI —7) 


By  subtracting  the  first  equation  of  (III-5)  from  (III-7),  we  see  that  the  equations  of  motion  of 
the  Earth-Moon  barycenter  relative  to  the  Sun  are 


,2  k 
d  xcs 

dt2 


/  M  \/m  x k  M  xk  \ 

=  _vM  L  +  _2  _j?  _es  _m  _ms 

yras  \  M  J  M  3  +  M  3  / 

s  \  c  r  „  c  r  / 
\  es  ms/ 


+  *k  +  ^{Fk  +  Fk} 

c 


k  =  1,  2,  3 


(III— 8) 
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where 


M.  /M  xk 


M  x.k  xk 


*  =  yMs  l  J 

.  .  s  \  c  r-  c  r.  r. 

3=1  \  je  jm  js, 


c  r.  r. 
jm  js 


,  k  =  1,  2,  3 


(III-9) 


Here  we  have  used  the  fact  that 


M  M  /  M  \ 

;  +  Ms  JjT  =  Ms  M~  [i  +  SC ') 
c  c  s 

M  M  /  M  \ 

+  1  s  M  s  M  \  M  /  ' 

r>  r>  >  c  ' 


Next,  to  obtain  the  equations  of  motion  of  the  Moon  about  the  Earth,  we  subtract  the  second 
equation  of  (III-5)  from  the  third  equation,  with  the  result  that 


d2xk  M  x k 

_ SI®  =  _yM  _c  _me 

.,2  7  s  M  3 

dt  s  r 

me 


,Bt**k,(r<-rFek)  ■  k  =  1’2'3  'UI-10> 

'  m  p  / 


where 


(k  x k 

¥  "HT 

r  r 
es  ms. 


.  N  M.  /x.k  xk 

*  =  yMs  l 

j  =  l  s  \jm  rje 


k  =  1,  2,  3  . 


(Ill-il) 


j  =  l  s  \rjm  rje/J 

Effects  contained  in  the  [(l/Mc)  (F^  +  F^)]  term  in  (III-8)  determined  in  this  report  are 

Rk  =  general  relativity  effect  [see  (IV-52)] 

S  =  second  harmonic  of  the  Sun  [see  (III— 50)  ]  . 

Effects  contained  in  the  [(l/M  )  Fk  -  {l/M  )  Fk]  term  in  (III- 10)  determined  in  this  report  are 

in  m  e  e 

k 

H  =  second  and  third  harmonics  of  the  Earth,  and  the  second  harmonic 
of  the  Moon. 

k  k 

If  we  let  F  and  F  denote  the  forces  acting  on  the  Earth  and  Moon  in  addition  to  those  enumev- 
e  m 

ated,  { III- 8 )  and  (III-10)  become 


d  X”  /  M  \  /  M  x  " 

~:f  - -™s  (“if) fef  + 


^  +  *k  Rk 

M  3  /  9  +  n 

c  r  / 
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M  '  e  m 
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k  =  1,2,  3 


(III-12) 


.2  k 
d  xme 


=  +  Bk  +  ^>k  +  Hk  +  Fk - —  Fk>) 

yi  s  M  3  M  1  e  / 

Q  r*  '  m  o  ' 
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If  we  desire  the  equations  of  motion  of  the  Earth-Moon  barycenter  to  have  the  appearance  of 
perturbed  elliptic  motion,  we  can  write  them  in  the  form 


d2xk 

- f  =  -yM 

dr 


/  M  x  x 
'  s'  r 

cs 


+  Ak  +  <I>k  +  Rk  +  SK 


+  j4-(F  k  +  Fk) 
M  '  e  m 
c 


1,2,3 


( III— 13) 


where 


,  /  M  \  /x k  M  xk  M 

A  -  /M  \1  +  M  H  3  m  3  M 

s '  \  r  „  c  r  _  c 
\  cs  es 


ms 


ms/ 


k  =  1,2,  3 


(III  — 14) 


C.  PLANETARY  PERTURBATIONS 

The  magnitude  of  the  acceleration  relative  to  the  Sun  that  the  Sun  gives  to  planet  p  is 


/  M  v  , 

yMs  1  +  M E)  4-  • 

'  s'  r 

ps 


(III— 1 5) 


The  magnitude  of  the  acceleration  relative  to  the  Sun  that  planet  j  gives  to  planet  p  at  its  closest 
approach  to  p  (assuming  that  the  Sun,  p  and  j  are  in  a  straight  line)  is 


M.  /  ,  ,  > 

yMs  m1  \~Z  ±_T 

S  V 


(III-16) 


where  the  plus  sign  ,s  used  if  j  is  between  p  and  the  Sun,  and  the  minus  sign  is  used  if  p  is  be¬ 
tween  j  and  the  Sun  [see  (III-2)  and  (III-3)).  Using  (III- 15)  and  (III-16),  the  equations  in  (IV-52) 

-4 

for  the  general  relativity  effect,  the  information  in  Table  I,  the  fact  that  (yM  )  =  2.96  X  10 
AU  /day  ,  and  the  following  discussion  of  the  Earth-Moon  dipolp  term,  we  derived  Table  II. 
According  to  (III-3),  we  can  write  the  effect  of  the  Earth  and  the  Moon  on  a  planet  as 

k 

ivi  ix  >r 

m  I  rr.  r» 
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M. 


a  =  yM  ~ 
em  s  M, 


M  /xk  xk 
_ e  /  ep  es 
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c  \  r  r 

mp  ms/ 


( III  — 17) 


The  effect  of  a  hypothetical  body  of  mass  Mc  =  Mg  +  Mm  at  the  Earth-Moon  barycenter  is 


nc  ■  ^Ms 


M  /Xk  ,k\ 

_ c  /  cp _ cs  \ 

M  l ~T  3  / 

s  \r  r„_/ 

\  cp  cs/ 


(III- 18) 


We  wish  to  determine  the  dipole  term 

Tk  =  a  k  -  q  k  . 

em  c 


(III-19) 


According  to  (III-6),  we  can  write 
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TABLE  I 

PLANETS  IN  THE  SOLAR  SYSTEM* 

(Note:  The  masses  of  the  outer  planets  include  those  of  their  satellites.) 


Mass 

Mean 

Period 

Inclination 

Planet 

(Sun  =  1) 

Distance  (AU) 

(years) 

Eccentricity 

(deg) 

1. 

Mercury 

1 .70  X  IQ'7 

0.387 

0.241 

0. 2056 

7.004 

2. 

Venus 

2.45X  10"6 

0.723 

0.615 

0.0068 

3.394 

3. 

Earth-Moon 

3.04X  10"6 

1.000 

1.000 

0.0167 

0 

barycenter 

4. 

Mars 

3.20  X 10~7 

1.524 

1.881 

0.0934 

1.850 

5. 

Jupiter 

9.55  X  10"4 

5.203 

1 1 . 862 

0.0484 

1.305 

6. 

Saturn 

2.85X  10-4 

9.539 

29.458 

0.0557 

2.490 

7. 

Uranus 

4.37  X  10'6 

19.18 

84.01 

0.0472 

0.773 

8. 

Neptune 

5. 18  X  10"5 

30.06 

164.79 

0.0086 

1.774 

9. 

Pluto 

2.78  X  10"6 

39.44 

247.69 

0.2502 

17. 170 

fSee  Ref. 5. 
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te)  (-6 


where 

3  £  £ 

y>  Xme  Xcs 

w  =  /, - 

s  r  r  „ 

^  me  cs 

with  exactly  similar  equations  holding  with  the  subscript  s  replaced  by  p 
we  then  have 
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where  the  Pf  are  Legendre  polynomials, 

[i/2]  i 

p,7,.i  V  (— l)l(2£-2i)i  „i-2i 

x  ^  “  ,|  i!  (£--i)!  (£  ~2i)' 

^  i=0 

The  first  few  Legendre  polynomials  are 

P0(Z)  =  1  ,  Pd(Z)  =  Z  ,  P2(Z)  =  |  Z2  -  j 


P3(Z)  =  f  Z3  -  |  Z 


P4(Z)=Tz4-T5z2+f  • 


From  (III— 2 1),  we  obtain 
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3  3 
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'W 
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where 


Q3f(Z)  = 


a1+a2+a3=£ 


P  (Z)  P  (Z)  Pa  (Z) 
“l  a2  a3 


According  to  ( III- 24),  we  see  that 


15  „2  3 


Q30(Z)  =  1  ,  Q31(Z)  =  3Z  ,  Q32(Z)  =  T  Z  -  f  • 


(Ill— 20) 

(III— 2 1) 

.  According  to  Ref.  6, 

(III— 22) 

(III— 23) 

(III— 24) 

(III— 25) 

(III— 26) 

(III— 27) 
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Finally,  in  formula  (III-19),  by  substituting  (III-25)  and  the  exactly  similar  equations  obtained 
by  replacing  the  subscript  s  by  p,  and  by  ignoring  powers  of  r  /r  and  r  /r  higher  than 

HiC  CS  1116  Cp 

the  second,  we  see  that 
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(111-28) 


The  last  column  in  Table  II  is  then  calculated  using  (III-28),  Table  I  and  the  following  ucu<t.‘ 
Mean  distance  of  Moon  from  Earth  =  384,400  km  =  0.0026  AU 


1 

82.31 


0.01215 


M  M 

jyp  =  1  -  -}yp  =  0.98785  .  (HI-29) 

c  c 

It  can  be  seen  that  the  effect  of  the  Earth  and  Moon  on  a  planet  can  be  written  as  (III- 18) 
rather  than  (III- 17),  even  in  the  case  of  Venus,  because  the  entry  in  Table  II  represents  the 
maximum  magnitude  of  the  dipole  term  when  the  planet  is  closest  to  t'ie  Earth,  and  this  maxi¬ 
mum  will  be  much  less  at  other  points  in  the  orbit.  Further,  the  sign  of  the  dipole  term  os¬ 
cillates  as  the  Moon  orbits  the  Earth. 

The  largest  satellite  of  Jupiter  is  Ganymede,  and  according  to  Ref.  7  we  have 

mass  of  Ganymede  _  1 

mass  of  Jupiter  "  12,300 

distance  from  Jupiter  =  7.156  X  10-3AU  .  (III-30) 

Thus,  (III-28)  implies  that  the  maximum  error  in  neglecting  the  displacement  of  Ganymede  from 
Jupiter  in  computing  the  perturbation  force  on  Mars  is 

j Tk|  =  1.45  X  10"17AU/day2 

which  is  three  orders  of  magnitude  less  than  the  maximum  effect  of  Pluto  on  Mars.  In  general, 
we  can  conclude  that  it  is  quite  accurate  to  assume  that  the  mass  of  an  outer  planet -satellite 
system  is  concentrated  at  the  center  of  mass  of  the  system  vrhen  determining  the  effect  of  the 
system  on  a  planet. 

The  total  mass  of  the  minor  planets  (asteroids)  is  estimated  to  be  about  3/lO,0OO  that  of  the 
8 

Earth.  It  would  be  rather  difficult  to  include  the  gravitational  effects  of  these  asteroids,  except 
for  some  of  the  larger  ones,  such  as  Ceres,  Pallas,  Juno  and  Vesta.  But  let  us  consider  the 
largest  asteroid,  Ceres.  According  to  Ref.  9,  we  have 

mass  of  Ceres  ,  ^  ,„-ll 

- -  r  o -  =  3.32  x  io 

mass  of  Sun 

mean  distance,  from  Sun  =  2.767  AU 
Thus,  by  (III- 16),  the  maximum  acceleration  of  Mars  due  to  Ceres  is 


a  =  5.0  X  10"15AU/day2 


26 


I 


L 


which  is  less  by  an  order  of  magnitude  than  the  effect  of  Pluto  on  Mars. 

The  maximum  acceleration  relative  to  the  Sun  that  a  distant  star  f  can  give  to  a  planet  p  of 
distance  A  from  the  Sun  is,  by  (III-16), 


fl  yms  <Mf/Ms>  [/rfs\2  / 

f  r2  !>V 

(— )2  +  ** 
L^fP/ 


yMs  (Mf/Mg) 


2yMg  (Mf/Ms)  A 


1  fs 


Here  we  have  assumed  that  r^g  =  r^.  +  A.  The  nearest  star  is  at  a  distance  of  4  light  years  = 

2.52  x  10^  AU  from  the  Sun.  Assuming  that  its  mass  is  the  same  as  the  Sun's,  we  have 

fif  =  3.70  X  10"20  AAU/day2 

where  A  ranges  between  0.387  ALT  for  Mercury  and  39.44  AU  for  Pluto.  To  obtain  the  accelera¬ 
tion  of  the  Moon  relative  to  the  Earth  due  to  the  star  f,  we  set  A  =  0.0026  AU.  Since  there  are  no 
really  massive  stars  in  the  neighborhood  of  the  Sun,  and  since  the  effect  of  a  star  on  the  accelera¬ 
tion  of  a  planet  relative  to  the  Sun  goes  down  as  the  cube  of  the  distance,  we  .an  feel  justified  in 
considering  the  solar  system  as  a  closed  system  in  discussing  the  orbital  motions  of  the  planets. 

The  effect  of  the  displacement  of  the  Moon  from  the  Earth  in  the  equations  of  motion  of  the 
Earth-Moon  barycentei  is  given  by  the  term  A  in  (III- 13) .  Inserting  (III-25)  in  formula  (III- 14) 

and  ignoring  powers  of  (r  /r  )  higher  than  the  second,  we  obtain 

me  cs 


Ak  =  yM 


/  M\MM/r  \2  . 

(.  c\  e  m(  me  )  1 

V  M  /  M  M  \  r  /  2 

s  c  c  cs  r  _ 
cs 


rxk  xk  .r 

-me  3w - —  (^w-2  -|) 

r  s  r  2  s  2 
me  cs 


A  simple  calculation  gives 

|  Ak|  <  2  X  10~10AU/day2 


(III— 3 1) 


(III— 32) 


so  that  it  is  important  to  retain  this  term  in  the  equations  of  motion  of  the  Earth-Moon  barycenter. 

Equations  (III-4)  for  the  motion  of  a  planet  in  the  case  of  a  planet  with  satellites  are  to  be  in¬ 
terpreted  as  the  equations  of  motion  of  the  center  of  mass  of  the  planet-satellite  system.  The 
error  in  these  equations  in  representing  the  motion  of  the  center  of  mass  is  given  by  a  term  simi- 
lar  to  the  term  A  of  (III-14).  Mercury  and  Venus  have  no  detectable  satellites,  and  the  satellites 
of  Mars  have  very  small  mass;  so  this  possible  error  is  only  of  concern  for  the  outer  planets. 
However,  by  (III-30)  and  (III-31),  the  error  in  ignoring  the  displacement  of  Ganymede  from  Jupi¬ 
ter  in  the  equations  of  motion  of  the  center  of  mass  of  the  Jovian  system  satisfies 

|Ak|  <  10'14AU/day2  . 

Further,  the  short  period  of  revolution  of  Ganymede  about  Jupiter  (7.15  days)  and  the  long  period 

k 

of  revolution  of  Jupiter  about  the  Sun  (11.86  years)  would  tend  to  cause  the  values  of  A  at  various 
times  io  cancel  each  other.  So,  in  general,  we  can  conclude  that  (III-4)  represents  the  motion  of 
the  center  of  mass  of  a  planet-satellite  system  in  the  case  of  an  outer  planet-satellite  system. 
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D.  SECOND  HARMONIC  OF  THE  GRAVITATIONAL  POTENTIAL  OF  THE  SUN 
12  3 

Let  (x  ,  x  ,  x  }  be  a  Cartesian  coordinate  system  with  origin  at  the  center  of  mass  of  a  body 
B  of  mass  M.  Let  dM  be  an  element  of  mass  of  the  body  B.  The  gravitational  potential  outside 
of  B  is  then 


U  =  -y  (Tf  dM ? ,f) 

•^B  [(x1  -  I1)2  +  (x2  -*)Z  +  (x3  -|3)2]4/ 


Let  r2  =  (x*)2  +  (x2)2  +  (x3)2.  We  have 

92U  ,  92U  92U  „ 

9(xV  9{x2)2  9(x3)2 


(III— 33) 


U  for  large  r 


(III— 3  4) 


We  introduce  spherical  coordinates  (r,  0,  (p)  by  the  formulas 


x  =  r  sin  9  cos  <p 


yi  =  r  sin  0  sin  cp)  0  <  r  <  «>,  0  ^  <p  <  Zn,  0  4  0  <  it 


(III— 35) 


x  =  r  cos  0 

Because  of  (III-34),  U  can  be  expanded  in  spherical  harmonics 


•vM  vM  V'  V'  ^>nh^COS 

u  =  -V  +  V  E  Z  <anhcoshp +  b  h  sinh?)  - - - - 


—  U  u  '  nh 

n=l  h=0 


(HI-36) 


where 


P  (Z)  =  P  (Z)  =  — —  (Z2  -  l)n  n  =  0,  1,  2 _ 

n0  n  2nn!  dZn 


dhP  (Z) 
n 


, ,  /7  a  r  { 

P  (Z)  =  (l-Z^r  - t- 

nn  dZn 


h  =  0,  1, .  . . , n 


(III— 3  7) 


(see  Refs.  10  and  11).  The  first  few  Legendre  polynomials  Pn(Z)  are  given  in  (III-24). 
Let  us  write 


r(x-|)  =  I  Yj  (**  -th: 


r=  l  <xj)2l 

Ljrl  J 

3  11 

f  =  E  u3>2 

•  i-4  ** 


(HI-38) 
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According  to  Ref.  6,  we  then  have 

co 

Tnr=a  ■  r  l  p,«i)  <?)' 

i= 0 

where  q  is  the  cosine  of  the  angle  between  the  vectors  r  and  f, 
3 

r  T 


(III— 3  9) 


<.=  z  4  4 

j=i 


(III-40) 


By  (III-33)  this  implies  that 

U  "  "T1  -  7  1  4  IijfB  p„«» 


n=l 


(IiI-41) 


Comparing  (III-41)  with  formula  (III-36),  we  see  that 
n 

M  Yi  (anh  coshl??  +  bnh  sinh<?>  pnh(c°S  0)  =  pn(q)  £ndM(£) 

h=0  "<B 

Since  the  origin  of  the  coordinate  system  is  at  the  center  of  mass  of  the  body,  we  have 

3 


(III-42) 


1TL  pi"i>  £dM(8  -Dr  m  *idM«  ’ 0  • 

3=1 


JB 


(III— 4  3) 


Thus,  the  summation  in  (III-36)  can  start  with  n  =  2, 

U  =  -iM  +  ^  J]  Yj  <anh  cosh'*  +  bnh  sinh<P) 


Pnh{coS  ^ 


r  ‘-J  ‘-i  '  nh 

n=2  h=0 


(III— 44) 


If  the  body  is  symmetric  about  a  line  through  its  center  of  mass,  and  we  choose  the  x-axis 
to  point  along  this  line,  (III-44)  reduces  to 


U  =  -XM  +  l!  y  ip  (cog  e) 
r  r  <->  n  n' 


(III— 45) 


n=2 


where  we  have  written  Jn  for  anQ.  In  the  case  of  the  Sun,  we  thus  suppose  that  the  gravitational 
potential  is  given  by 

yM 


yM  /Re\  2/S,\2 

U  - - -  +  - - 

r  r  , 


ft)2©1  [!(#-! 


(111-46) 


where  R  =  6.96  X  10s  km  is  the  radius  of  the  Sun,  and  the  coefficient  S7  of  the  second  harmonic 

of  the  Sun  is  to  be  determined  by  its  effect  on  the  motion  of  a  planet.  We  may  assume  that  the 

x3-axis  of  symmetry  of  the  Sun  points  along  the  axis  of  rotation  of  the  Sun. 

12  3 

Let  (X  ,  X  ,  X  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Sun  in  which 
the  equations  of  motion  are  expressed.  We  have 
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We  recall  that  S£  is  the  second  harmonic  of  the  Sun's  gravitational  potential  and  has  the  dimen¬ 
sions  of  a  length  squared,  and  that  Rg  is  the  equatorial  radius  of  the  Sun.  The  quantities 
C3i  (i  =  1,  2,  3)  are  determined  in  Appendix  C. 

E.  HIGHER  HARMONICS  IN  THE  GRAVITATIONAL  POTENTIALS 
OF  THE  EARTH  AND  MOON 


The  purpose  of  this  section  is  to  derive  the  force  on  the  Moon  due  to  the  Earth,  .considering 
terms  up  to  the  third  harmonic  in  the  Earth's  gravitational  potential  and  up  to  the  second  harmonic 
in  the  Moon's  gravitational  potential. 

We  shall  assume  that  the  Earth  is  symmetric  about  its  axis  of  rotation.  Let  the  coordinate 
12  3 

system  (x  ,  x  ,  x  ),  with  origin  at  the  center  of  mass  of  the  Earth,  be  referred  to  the  true  equinox 

3 

and  equator  of  date,  so  that  the  x  -axis  points  along  this  axis  of  rotation.  Then,  by  (III-45)  and 
(III-24),  we  have 


30 


{Ill— 52) 


yM  yM  J-)  ^  ?  n  yM  e  3 

u  =  -  +  — %r-  (4  cos2  0  -  b  +  — (f  '’os3  0  -  |  cos  ©) 


where,  by  (III-35), 


cos  0  =  — 
r 


According  to  Ref.  12,  we  have 
J 


(III— 53) 


R 


^  =  1.0827  X  10"3 


V  o  / 

—A  =  _2.4  X  10“° 

R- 

e 


(III— 54) 


where  R  =  6378.17  km  is  the  equatorial  radius  of  the  Earth. 

6  1  o  o 

Let  (X  ,  X  ,  X  )  be  the  coordinate  system  referred  to  the  mean  equinox  and  equator  of 
1950.0,  the  reference  system  in  which  we  are  integrating  the  equations  of  motion.  Choc  sing 
the  origin  of  this  coordinate  system  to  be  at  the  center  of  mass  of  the  Earth,  we  can  write 


=  Z  X< 


t=  1 


xi  ■  Z  A« 


f=l 


j  =  1,2,3 


(III-55) 


where  the  orthogonal  matrix  (A^)  is  given  in  Appendix  A.  Then,  by  (111-52)  through  (HI-55)  and 
(III-48),  the  components  of  force  (F1,  F2,  F3)  on  a  particle -of  mass  M  due  to  the  Earth  in  the 


12  3 

coordinate  system  (X  ,  X  ,  X  )  are 

rk 


Fk(X)  =  - 


[xk  ,15 


r 

yM  MJ, 


(y  COS2  0  -  |)  -  3A3k  cos  0 


el 


£  (25  cos3 e  -  cose) 

r  2  2 


-A3k  (i§  cos2e  -  |) 


k  =  1,.2,  3 


(III-56) 


Suppose  that  B  is  an  extended  physical  body  of  mass  M  (the  Moon  in  our  case),  and  let 
(X1,  X2,  X3)  be  the  coordinates  of  its  center  of  mass.  Let  (51,  |2,  |3)  be  the  coordinates  of  a  mass 
element  dM  in  B  relative  to  the  center  of  mass  of  B.  By  (III-56),  the  force  Fg  on  B  due  to  the 
gravitational  field  of  the  Earth  has  components 


Fk(X)  = 


B 


Fk(X  +  |)  =  Fk(X)  +  Ck 


(III— 5  7) 


where 
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We  write 


where 


where 


Ck  =  [Fk<X  +  Q  -  Fk(X)] 


Ck  =  C  k  +  C  k  +  C  k 


Cf  -  ~vMe 


k  . 


err  [x_t 

lr(X  + 


X 


r{X  +  S)3  r(X)3 


dM{?) 


C  k  =  yMeJ2  [Gk{X  +  8-  Gk(X)]  dM($ 

C  k  =  yMeJ3  [Hk(X  +  |)  -  Hk(X)]  dM{0 


r  3 


r(X  -  «)  = 


£  (xj  +  ij)2 

Lj=l 


1/2 


Xk  ,15 _ 2  -  3.  3A3kCOS0 


(X)  =  (y  cos  0  -  |)  - 


Hk(X)  =  -^r  (¥  cos3  ©  -  y  cos  e)  — cosZ  9  ~  f) 
r  r 


Equation  (III-39)  implies  that 

r(X+T)  =  r  £  (_1)  PC(q)  ^ 


£=0 


(III— 58) 


(III— 59) 


(III— 60) 


(III-61) 


r  ~ 


£  = 


where  the  P^(q)  are  Legendre  polynomials  and  where 

3  1 1/2 

£  (Xj)2 
Lj  =  l 

r  3  ,1/2 

£  (^)2I 

Lj=l 


q  = 


y  x^  ll 

L>  r  £ 

j=i 


(III-62) 


Since  P.(q)  only  contains  even  powers  of  q  for  £  even,  and  odd  powers  of  q  for  t  odd,  expres- 

*  12  3 

sion  (III-61)  contains  no  square  roots  of  the  quantities  (£  ,  £  ,  £  ),  only  products  and  powers.  In 

(  12  3 

fact,  P^(q)  f  is  a  homogeneous  polynomial  in  (|  ,  £  ,  %  )  of  degree  £  with  coefficients  depending 
on  (xVr,  X2/r,  X3/r).  Equation  (111-61)  implies  that 
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1  _  _1_ 
r(X  +  |)n  r11 


2 

£=G 


Qne(q>  =  E  pa  (q)---pa  <q>  (HI-63) 

a.+. . .  +a  =£  1  n 

1  n 

(  12  3 

where  Qn((q)f  is  a  homogeneous  polynomial  in  (|  ,  |  |  )  of  degree  (  with  coefficients  depending 

on  (X1/ r,  X2/r,  X3/r). 

If  we  insert  the  above  expression  for  l/r(X  +  £  )n  into  the  integrands  in  (III-59),  we  obtain 


C 


k 

1 


yM 
_ _ e 

2 

r 


R^U.  f)  dM(«) 


R^U,  £)  dM(|) 


rMeJ3 


(III— 64) 


where  the  R^  [£,  (X/r)]  are  homogeneous  polynomials  in  (I1,  ,  |3)  of  degree  l  with  coefficients 

depending  on  (xVr,  X2/r,  X3/r).  Actually,  the  above  series  start  with  (  =  1,  because  the 

12  3 

integrands  in  (III-59)  are  the  differences  of  functions  evaluated  at  (X  +  £)  and  X.  Since  (4  ,  4  ,  4  ) 
(0,  0,  0)  is  the  center  of  mass  of  the  body, 


j  =  1.  2,  3 


(III— 65) 


so  we  may  assume  that  the  series  ( III- 64)  start  with  1  =  2.  The  integrals  ///g  Rjk  dM  (4)  involve 
the  moments  and  products  of  inertia  of  the  body,  while  the  integrals  fffg  R^dM(4)  ( l  >  2)  depend 
on  the  higher  moments  of  the  body.  We  ignore  these  higher  moments  which,  by  the  discussion 
in  Ref.  13,  is  equivalent  to  ignoring  harmonics  higher  than  the  second  in  the  body's  (=  Moon's) 
gravitational  potential.  We  can  therefore  assume  that 


rk_ 

S  "  4 

r 


B 


R^d,  f)  dMd) 


yMeJ2 


yMeJ3 


7 

r 


wB 


R2k2U. 

R32^' 


7>  dM(|) 
f)  dM(4) 


(III-66) 


1c 

The  effect  of  the  third  harmonic  of  the  Earth's  gravitational  potential  in  the  force  F  (X)  of 

(III-56)  is  of  the  order  J,/r  .  Thus,  to  the  accuracy  to  which  we  are  working  in  this  section, 

k  3  k  /  6  /  7 

we  can  assume  that  C2  and  are  zero,  because  the  coefficients  multiplying  1/r  and  1/r  in 

C  k  and  C  k  will  involve  J.,  and  times  similar  constants  associated  with  the  second  harmonic 
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-4; 

»  i 

\  j 


) 


i 


lc  k 

of  the  Moon's  gravitational  potential.  {The  fact  that  C ^  and  are  no  larger  than  the  effect  of 
the  fourth  harmonic  of  the  Earth  was  actually  checked  by  direct  computation,  but  the  calcula¬ 
tions  are  too  lengthy  to  reproduce  in  this  report. )  Finally,  examining  expression  (111-59)  for 

k  k  . 

,  and  noting  (III-33)  and  (III-48),  we  may  assume  that  is  minus  the  force  due  to  the  second 

harmonic  of  the  Moon's  gravitational  potential  acting  on  the  Earth  (where  the  Earth  is  considered 

to  be  a  point  mass). 
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The  Moon  is  approximately  an  ellipsoid  with  three  unequal  axes.  Let  (x  ,  x  ,  x  )  be  a  co- 

3 

ordinate  system  with  origin  at  the  center  of  mass  of  the  Moon  such  that  the  x  -axis  points  to- 

ward  the  north  pole  of  the  Moon  (which  is  one  of  the  axes  of  the  ellipsoid),  the  x  -axis  points 

2 

along  the  axis  of  the  ellipsoid  pointed  in  the  direction  of  the  Earth,  and  the  x  -axis  completes 

the  right-hand  system.  Let  L  be  the  moment  of  inertia  of  the  Moon  with  respect  to  the  x^-axis. 

J  4  2  3 

We  may  assume  that  the  products  of  inertia  with  respect  to  the  (x  ,  x  ,  x  )  frame  are  zero,  so 
that 


JB 


E  <*£>2 


dMU)  =  L 


j  =  1,  2,  3 


dM(l)  =  °  ’  i^j 


(III- 67 ) 
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Let  a,  b,  c  be  the  axes  of  the  ellipsoidal  Moon  in  the  x  ,  x  ,  x  directions.  According  to  Ref.  14, 
we  have 


b  +  c 


=  1737.9  km 


a  -  c  =  1.09  km 
a  —  b  =  0.36  km 
I, 


(III— 6  8) 


bZM 


=  0.397 


m 


h~lZ 


=  0.000420 


=  0.000628 


h-h 


=  0.000208 


(III-69) 


The  above  values  of  the  moments  of  inertia  of  the  Moon  (wnich  determine  the  second  harmonic 
of  the  Moon's  gravitational  potential)  were  obtained  from  the  observed  shape  and  physical  libra- 
tion  of  the  Moon.  It  is  to  be  expected  that  in  the  near  future  the  second  and  higher  harmonics 
of  the  Moon's  gravitational  potential  will  be  accurately  determined  by  placing  an  artificial  satel¬ 
lite  in  orbit  about  the  Moon. 

According  to  (III-24),  (III-42)  and  (III-44),  the  second  harmonic  of  the  Moon's  gravitational 
potential  is 
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u. 


=  "5  ffl 


P7(q)  rdMU) 


y 

3 

r 


i»  j 


Equations  (III-67)  give 
3 


"ryS  h3^)2]!!!  (4S)2dMlil 


2r 

_  _y 

4r 


j=i 

3 


2  [‘-’(rVlfE  ■»->,)  • 

3  =  1  'f7-j  7 


<3  *  2 

Using  the  fact  that  2  (xVr)  =  1,  this  can  be  put  in  the  form 
3=1 


“2-Ji  2  [!  (t)2  -  I 


3  =  1 
3 


I. 

3 


_y 

3 


2  Kt) 


3  =  2 


3\2 


(L-Il) 


Introducing  polar  coordinates  (III-35),  the  first  line  of  (III- 70)  implies  that 

U2  =  (|  (1^  cos2  ©  +  i,  sin2  <p)  sin2  0  +  -|  I3  cos2  0  -  |(I^  +  I2  +  I3)) 

r 

=  -3  {^3  -  I  (I1  +  I2)1  f|  cos2e  -  |l  +  3  di  -y  cos  2^  Sin2©} 
r 

Comparing  this  expression  with  (III-37),  (III-41)  ynd  (III-42),  we  see  that 
yM 


U->  =  — [a2QP20(cos  0)  +  a22  cos  2<p  P22(cos  0)  J 
r 


where 


M  a.-  =  I_  —  =•  (I,  +  I9) 
m  20  3  2  1  2' 


(III— 7  0) 


(III— 7  1) 


(III-72) 


M  a„  =  7  (I,  -  I,) 
m  22  4  1  2' 


which  implies  that 


IQ  -  I .  =  4M  a,, 
2  1  m  22 


I,  -  I .  =  M  (a9r,  —  2a99) 
3  1  m  20  22' 


(III— 73) 


(III  —  7  4) 


We  shall  use  expression  (III-70)  for  the  second  harmonic  of  the  Moon's  gravitational  potential. 
We  have  derived  (III-72)  through  ( III- 74)  so  that  we  can  determine  improved  values  of  the 
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quantities  (I^  - 1^)  and  (I^  —  1^)  in  (III-70)  from  a  possible  future  publication  of  improved  values 
of  a20  and  a22. 

Let  (X*,  X2,  X3)  be  the  coordinates  of  the  center  of  mass  of  the  Moon  relative  to  the  center 

of  mass  of  the  Earth  in  the  coordinate  system  in  which  we  are  integrating  the  equations  of  motion. 
12  3 

Let  (x  ,  x  ,  x  )  be  the  coordinates  of  the  center  of  mass  of  the  Earth  relative  to  the  center  of 
mass  of  the  Moon  in  the  coordinate  system  used  in  formula  (III-70).  Then,  the  relation  between 
(X1,  X2,  X2)  and  (x*,  x2,  x3)  is  given  by 


Bj(x' 


£=1 


xj  *  -  Z  B(]*' 


£=1 


j  =  1,  2,  3 


(III— 7  5) 


where  the  orthogonal  matrix  B. .  is  determined  in  Appendix  B. 
k  *3  12  3 

Let  be  the  components  in  the  (X  ,  X  ,  X  )  coordinate  system  of  minus  the  force  due  to 
the  second  harmonic  of  the  Moon's  gravitational  potential  acting  on  the  Earth  (where  the  Earth 
is  considered  to  be  a  point  mass).  Then,  by  (III-48),  (III-70)  and  (II1-75),  we  have 

au. 


Ck  =  -M 
1  e 


yM 

r 


ax- 

3 


I1  Z  W--L)  [“ (^  D2  -  |) -3D.B 
4  j  1'  [  r  2  3  2'  j  jk 

j=2 


(III— 7  6) 


where 


D.  =- 
3 


2d  -  y  b  — 

r  “  Lj  j£  r 


£=1 


(III— 77) 


Formulas  (III-56),  (III-57),  (III-58)  and  (III-76)  combine  to  give  that  the  force  on  the  Moon  due  to 
the  Earth  is 


,  yM  M  Xk  yM  M  /R  \2  J?  rvk  ,  _ 

„k  _ e  m  ,  '  e  m  /  e\  2  X  ,15  2  _  3.  _A  _ 

Fm  r3  +  r2  \  r  }  R2  [  r  (  2  COS  0  2}  3A3k  C0S  0 


?MeMm  f^n\Z 


v  Ii“I 
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1  |X-(TD32-|-3DjB3k 
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e  m 


/Re\3  J3  fxk  ,35  3  -  15 

^3  [—  (T  COS  0  "  T  COS0) 


IS  2  S 

-A3k(T  cos  9  -  V 


k  =  1,  2,  3 


(III-78) 


where  R  and  Rm  are  the  radii  of  the  Earth  and  Moon,  and  where  we  have  assumed  that  C,  =  0, 

i.  e  k  ^ 

Cj  =  0-  The  force  F  on  the  Earth  due  to  the  Moon  is  minus  the  force  on  the  Moon  due  to  the 
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Earth,  Fk  =  —  F^.  Equations  (III- 10)  for  the  motion  of  the  Moon  relative  to  the  Earth  were 
derived  by  subtracting  the  equations  of  motion  of  the  Earth  from  the  equations  of  motion  of  the 
Moon.  Thus,  the  Hk  term  in  (III-12)  is 


Hk  = 
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where  we  have  used  the  notation  of  Sec.  III-B  with  x 


me 


r  =  r,  and  M  =  M  +  M 
me  c  e  m 


(III— 79) 
The 


constants  (Re,  J2,  J3)  and  (Rm,  I1, 12, 13)  are  given  in  (111-54),  (III-68)  and  (III-69).  (We  must 
evidently  assume  that  Rm  =  b.)  The  matrices  (A^)  and  (B.j)  are  determined  in  Appendices  A 
and  B.  By  (111-53)  and  (III-77),  the  quantities  cos  0  and  in  (111-79)  are 
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E».  =  T.  B. 
3  u  3 


3  xf 
‘  me 


£=1 


j£  r 


me 


3  =  2,3 


(III— 80) 


Table  III  is  constructed  using  the  expressions  for  the  forces  perturbing  the  motion  of  the 
Moon  relative  to  the  Earth  given  in  (III-ll),  (III-12)  and  (III-79). 


TABLE  III 

MAXIMUM  ACCELERATIONS  OF  THE  MOON  RELATIVE  TO  THE  EARTH  (AU/DAY2) 
(Note:  A  constant  acceleration  of  10-^  AU/day^  will  move  a  body  1  km  in  10  years.) 

Due  to 

Acceleration 

Due  to 

Acceleration 

Earth 

1 .33  X  10"4 

Mars 

3.45X10"12 

Sun 

1.55X  10"6 

Jupiter 

1 . 98  X  1 0" 1 1 

2nd  harmonic  of  Earth 

-10 

2.40X10 

Saturn 

7.C5X  10"13 

2nd  harmonic  of  Moon 

3.00X-.0"12 

Uranus 

1 . 1 2  X  1 0** 1 4 

3rd  harmonic  of  Earth 

2.00X10'14 

Neptune 

3.25X10"15 

Mercury 

-12 

1.14X10 

Pluto 

7.53  X  10“17 

Venus 

I.80X  10“10 
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The  higher  harmonics  of  the  Earth's  gravitational  potential  are  known  quite  accurately  be¬ 
cause  of  their  effects  on  the  motions  of  artificial  Earth  satellites.  If  the  higher  harmonics  of 
the  Moon's  gravitational  potential  are  similarly  determined  by  placing  an  artificial  satellite  in 
orbit  about  the  Moon,  the  effect  of  both  the  second  and  third  harmonics  of  the  Moon  can  be  m- 
eluded  in  the  H  term  above.  In  this  case,  the  second,  third  and  fourth  harmonics  of  the  Earth, 
the  interaction  between  the  second  harmonics  of  the  Earth  and  Moon  [the  C  ^  term  in  (111-66)], 
and  the  effect  of  the  terms  in  the  gravitational  potential  of  the  Earth  which  arise  from  asym¬ 
metries  about  the  north-south  axis  (tesseral  harmonics)  can  all  be  included. 

The  effect  of  tidal  friction  on  the  motion  of  the  Moon  is  small,  since  it  is  estimated  that  the 

-4  15 

increase  in  the  sidereal  day  as  a  result  of  tidal  action  is  7.2  X  10  sec  per  century.  However, 
the  effect  of  tidal  friction  should  be  included  in  the  equations  of  motion  of  the  Moon,  if  harmonics 
in  ihe  gravitational  potentials  of  the  Earth  and  Moon  higher  than  those  considered  in  this  report 
are  included  in  the  equations  of  motion. 
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IV.  GENERAL  RELATIVISTIC  EFFECT 


A.  MATHEMATICAL  FORMALISM  OF  EINSTEIN'S  GENERAL  THEORY 
OF  RELATIVITY 

Euclidean  n-space  Kn  is  defined  as  the  set  of  all  n-tuples  (x*.  ....  x11)  of  real  numbers 
X"*  (j  =  1,  . . . ,  n).  An  open  ball  of  radius  r  about  a  point  xq  =  (xj ,  ....  x” )  in  !Rn  is  the  set  of  all 

1  /~n  ‘  •  2 

points  x  =  (x1,  ....  xn)  in  such  that  /  2  (x^  —  x**)  <  r.  A  manifold  Mn  of  dimension  n  is  a 

v  j---l  0 

J  n 

separable  connected  Hausdorff  space  such  ‘hat  each  point  in  M  has  a  neighborhood  which  is 

homeomorphic  to  an  open  ball  in  $1?.^  That  is,  each  point  p  in  Mn  has  a  neighborhood  U  in  Mn 

such  that  there  is  a  one-to-one  bicontinuous  map  f:U  —  B  onto  an  open  ball  B  in  £Rn.  The  map 

f  puts  a  coordinate  system  on  U  in  the  sense  that  io  each  point  q  in  U  is  associated  the  co- 
1  n 

ordinates  (x  ,  .  .  . ,  x  )  =  f(q).  Suppose  that  f^U^  -*  and  f2:U2  B2  are  two  coordinate  systems 
such  that  the  intersection  of  U .  and  U-  is  not  empty.  Then  the  change  of  coordinates  is  given  by 

*  L 

the  map  f2  °  :B^  —  B^,  as  sketched  in  Fig  2.  A  differentiable  manifold  of  dimension  n  is  a 

manifold  of  dimension  n  which  has  a  covering  by  coordinate  neighborhoods  such  that  the  co¬ 
ordinate  transformations  are  infinitely  differentiable.  The  manifold  is  analytic  if  the  coordinate 
transformations  are  analytic. 


The  simplest  example  of  an  n-dimensiona\  manifold  is  Euclidean  n-space  itself.  Examples 
of  two  dimensional  manifolds  are  provided  by  surfaces  in  Euclidean  three-space,  such  as  the 
cylinder,  torus  and  sphere.  A  manifold  can  be  defined  without  any  reference  to  a  higher  dimen¬ 
sional  Euclidean  space.  Roughly  speaking,  one  might  imagine  that  a  manifold  is  a  space  which 
can  curve  back  on  itself  in  the  large,  tut  which  locally  looks  like  Euclidean  space. 

Let  (x*,  .  .  .  ,  x11)  and  (y*,  .  . . ,  yn)  be  two  overlappmg  coordinate  systems  on  a  differentiable 
manifold  Mn  of  dimension  n.  A  tensor  T  contravariant  of  order  p  and  covariant  of  order  q  is 
expressed  in  these  two  coordinate  systems  in  the  form^ 

tSee  Ref.  16  for  definitions  of  these  topological  concepts. 

t  See  Ref.  17  for  a  rigorous  abstract  definition  of  tensors  on  a  differerdiable  manifold. 
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p....p  t».  v  a  n 

T  =  T  1  /dx  4®...®dx  q®-~ 

r"  q  &x  1  avP 
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Lr>  ^  1  &  n  a 

T  =  S„  *  „  p  dy  1  ®. . .  ®  dy  q  ®  — - 
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(IV  -1) 


9y 


9y 


Here,  we  use  the  Einstein  summation  convention  in  which  x-epeated  upper  and  lower  indices  are 
summed.  The  components  of  the  tensor  T  transform  according  to  the  tensor  rule  of  trans¬ 
formation: 

v,  v  a.  a 


a 


.  .  ap  _  T  ‘  ‘  pp  3x  1  3x  q  9y  1  ay  p 

V-*a  - K  ~~~ 

1  q  q  9y  ay  q 


o  x 


ax^P 


(IV-2) 


In  the  general  theory  of  relativity,  the  space-time  universe  is  imagined  to  be  a  four¬ 
dimensional  differentiable  {or  perhaps  analytic)  manifold.  The  gravitational  potential  in  the 

2 

space-time  universe  is  given  by  a  symmetric  hyperbolic  covariant  tensor  of  order  two  ds  = 
g^dx*1  ®  dx^,  called  the  metric  tensor.  Symmetric  means  that  g^^  =  g,^,  and  hyperbolic  means 
that  for  each  point  of  space-time  universe  there  is  a  coordinate  system  (x®,  x1,  x^,  x3)  such  that 

ds2  =  dx°  <g)  dx°  -  dx1  ®  dx1  -  dx2  ®  dx2  -  dx3  (x)  dx3  (IV-3) 


at  that  point. 

Let  (ga^)  denote  the  matrix  inverse  to  the  matrix  (g  ).  The  Christoffel  symbols  are  then 
defined  by 


ra  _  1  af}(dgp\i  ,  9g0r  _  8gpi»\ 

V  "  2  g  'ax'"  zJ  I 


3xh 


a/ 


(IV -4) 


(see  Ref.  18).  If  we  think  of  the  metric  tensor  g  as  being  the  gravitational  potential,  then  we 

[Xls 

should  think  of  the  Christoffel  symbols  as  being  the  gravitational  field.  The  Riemann  curvature 
tensor  is  defined  by 


R 


a 


ax" 


OF 


I  p  ^  p  ® 

p/3  A  V  \IV  \($ 


(IV-5) 


We  further  define 


R  =  B.P  . 
pr  pr/3 

R  =  g^R^  (IV-6) 

(see  Ref.  19). 

It  is  then  postulated  that  the  gravitational  potential  in  the  space-time  universe  satisfies  the 
Einstein  field  equations 

R  -4-g  R  =  -  «T  (IV-7) 

pr  2  6pr  pr  '  ' 

where  *  is  a  constant  and  where  T  is  a  tensor  defined  in  terms  of  the  distribution  of  matter 
and  energy  in  the  space -time  universe.  Multiplying  both  sides  of  this  equation  by  g1^ ,  and 
summing  on  p  and  v,  we  see  that 
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R  =  k  T  .  (IV -8) 

If  there  is  no  matter  at  a  given  point  of  the  space-time  universe,  then  T  -  0  and  by  (IV-7)  and 

fit' 

(IV -8)  the  Einstein  field  equations  become 

V  =  0  <IV-9) 

at  this  point. 

The  above  equations  involving  the  metric,  Christoffel  symbols  and  Riemannian  curvature 
can  be  expressed  in  more  abstract  differential-geometric  terms.  This  abstract  approach  would 
be  appropriate  for  a  discussion  of  the  space -time  universe  in  the  large.  For  local  discussions, 
the  formulation  in  terms  of  local  coordinate  systems  is  sufficient. 

A  curve,  or  world  line,  in  the  space-time  universe  M  is  a  map  y:  (a,  b ]  "  M  of  an  interval 
[a,  b]  in  the  real  numbers  into  M  (see  Fig.  3).  In  a  coordinate  system  (x^,  x*,  x2,  x3)  on  a  co¬ 
ordinate  neighborhood  U  in  M,  the  curve  can  be  written  in  the  form 

x(i  =  yH(s)  ,  s  e  (a,  b]  .  (IV -10) 

The  tangent  vector  to  the  curve  is  then 


(IV-11) 


A  vector  X  =  X^O/ax11)  at  a  point  is  said  to  be  time-like  if  g  ,  X^X^  >  0,  null  if  g  X^X17  =  0,  and 

H  1/ 

space-like  if  g  X'x  <  0.  The  path  of  a  light  ray  through  the  space-time  universe  has  null 

flK 

tangent  vector,  while  the  path  of  a  material  body  through  the  space-time  univeise  has  time-like 
tangent  vector.  A  curve  in  the  space-time  universe  with  space-like  tangent  vector  has  no  physi¬ 
cal  interpretation. 

An  observer  in  the  space-time  universe  follows  a  time-like  world  line  y  through  the  space- 
time  universe.  Suppose  that  this  observer  possesses  an  atomic  clock  and  that  he  defines  a  sec¬ 
ond  of  time  to  be  a  certain  number  of  oscillations  of  this  atomic  clock.  Then,  in  traveling  along 
Irs  world  line  y  from  the  point  y(a)  to  the  point  y(b)  in  the  space -time  universe.,  it  is  postulated 
that  the  observer  will  see  that  the  number  of  elapsed  seconds  is  given  by  the  proper  time  integral 


Fig. 3.  Curve  on  a  manifold. 
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(IV-12) 


_  1  fb 

Tab  '  5  4 

where  c  is  a  constant,  dependent  upon  the  specific  chemical  element  whose  atomic  oscillations 
run  the  atomic  clock  and  upon  the  number  of  oscillations  defined  to  be  in  a  second.  The  above 
integral  depends  only  on  the  world  line  y  and  the  gravitational  potential  g  .  The  world  line  de- 

(XV 

pends  on  whether  the  observer  is  accelerated,  etc.  Thus,  it  is  only  reasonable  to  postulate 
(IV-12)  for  the  rate  of  a  clock  for  an  ideal  atomic  clock,  since  the  effect  of  impulses  on  the  rate 
of  a  mechanical  clock  would  depend  on  the  details  of  its  construction. 

Let  y:[a,  b]  —  M  be  a  null  or  time-like  curve  in  the  space-time  universe  M.  The  length  of 
y  is  defined  as 


[  dy^  dy^ 
ds  ds 


ds 


L(y)  = 


dx*1  dx^ 
ds  ds 


ds 


(IV -13) 


The  curve  y  is  a  geodesic  if  L(y)  is  a  minimum  for  all  nearby  curves  joining  y(a)  and  y(b)  in  M. 
If  the  parameter  s  satisfies  g  (dx^/ds)  (dxVds)  -  constant  along  the  curve,  then  a  null  or  time- 
like  geodesic  satisfies  the  differential  equations 


d2x^  r/3  dx^  dx^  _  Q 
^2  pr  ds  ds 


(IV-14) 


(see  Ref.  20).  It  is  postulated  that  the  path  followed  by  a  particle  of  negligible  mass  through  the 
space-time  universe,  subject  to  no  force  except  that  due  to  the  gravitational  potential  in  the 
space-time  universe,  is  a  time-like  geodesic. 

In  order  to  employ  the  theoretical  facade  outlined  above  ir.  concrete  situations,  we  make  the 

12  3 

following  comments.  First,  in  an  inertial  coordinate  system  (t,  x  ,x  ,x  )  of  special  relativity 
far  removed  from  ponderable  matter,  the  gravitational  potential  should  assume  the  form 


ds2  =  c2dt2  -  (dx1)2  -  (dx2)2  -  (dx3)2  .  (IV-15) 

2 

Here,  c  is  the  velocity  of  light,  and  dt  is  short  for  dt  <g)  dt.  In  the  spherical  coordinate  system 
(t,  r,  0,  (p ),  defined  by 


t  =  t 
1 

x  =  r  sin  6  cos  <p 
2 

x  =  r  sin  0  sin  <p 


x 


3 


=  r  cos  © 


(IV-16) 


the  metric  (IV-14)  becomes 

,  2  2, ,2  ,  2  2  ,02  2  .  2_  ,  2  .... 

ds  =  c  dt  -  dr  -  r  d©  -  r  sin  ©  dip  ,  (IV-17) 

Second,  in  a  general  relativistic  coordinate  system  which  closely  approximates  a  Newtonian 

inertial  coordinate  system,  the  Newtonian  expression  for  the  motion  in  a  weak  gravitational 

potential  U  of  a  particle  of  small  mass  with  velocity  small  relative  to  the  velocity  of  light  should 

be  approximately  the  same  as  the  general  relativistic  expression.  We  therefore  have  approxi- 
22 

mately 
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ds2  =  (c2  4  2U)  dt2  -  (dx1)2  -  (dx2)2  -  (dx3)2  .  (IV-18) 

Here,  we  assume  that  U  goes  to  zero  at  spatial  infinity  and  that  the  Newtonian  acceleration  of 

a  particle  is  —grad  U.  {The  convention  in  Ref.  22  is  that  the  Newtonian  acceleration  of  a  particle 

is  t-grad  U.)  Using  the  fact  that  (IV-18)  satisfies  the  Einstein  field  equations,  a  more  exact  ex- 

23 

pression  for  the  metric  is 

ds2  =  (c2  +  2U)  dt2  -  (l  -  -|r)  [(dx1)2  +  (dx2)2  +  (dx3)2]  .  (IV-19) 

'  c  ^ 


B.  MOTION  OF  A  PLANET  OF  SMALL  MASS  IN  THE  GRAVITATIONAL 
FIELD  OF  THE  SUN 


Let  V  be  a  coordinate  neighborhood  in  the  space -time  universe  containing  an  isolated 

12  3 

spherically  symmetric  body,  and  suppose  that  (t,  x  ,  x  ,  x  )  are  coordinates  on  V  such  that  the 
center  of  the  body  follows  the  world  line  given  by  x-*  =  0  (j  -  1,  2,  3).  The  gravitational  potential 
on  V  may  be  written 

ds2  =  gOQdt  ®  dt  +  g  .dt  ®  dx1  +  g^dx1  <x)  dt  +  g^dx1  ®  dx-*  .  (IV-20) 

Here,  and  in  the  following,  we  assume  that  Roman  indices  i,  j, . . .  take  on  only  the  values  1,  2,  3, 
Since  the  body  which  generates  the  gravitational  potential  (IV-20)  is  spherically  symmetric  and 
isolated,  we  may  suppose  that 

(1)  the  gravitational  potential  is  static,  i.e.,  the  components  of  the  metric 
tensor  do  not  depend  on  the  variable  t; 

(2)  the  line  element  (IV-20)  does  not  change  its  form  under  a  rotation  of  the 
coordinate  axes  (x1,  x2,  x3); 

(3)  at  a  large  spatial  distance  r  =  (x1)2  +  (x2)2  +  (x3)2  from  the  body,  (IV-20) 

approaches  the  value  (IV- 15). 


From  these  assumptions  and  the  fact  that  outside  the  body  the  Einstein  field  equations  (IV-9)  are 
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satisfied,  it  follows  that  there  is  a  coordinate  system  (L.,  x...,  x^,  x*)  on  V  such  that  outside 
the  body  the  metric  tensor  has  components24 
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(IV -21) 
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Here,  the  Kronecker  delta  6.^  is  defined  by  (11-60),  r...  =  ^(x*)  +  (x,.t)  +  (xi)  ,  c  is  the  velocity 
of  light  at  a  large  spatial  distance  from  the  body,  and  a  is  a  constant.  Comparison  of  (IV-21) 
and  (IV-18)  with  U  =  -(yM/r...)  shows  that 


a 


(IV -22) 


where  y  is  the  gravitational  constant,  and  M  is  the  mass  of  the  body.  The  constant  a  has  the 
dimensions  of  a  length  and  is  much  smaller  than  the  geometric  radius  of  the  body  (in  the  case  of 

the  Sun,  a  -  1.48  km).  In  the  spherical  coordinate  system  (t,.,,  r...,  0..,,  (pif)  defined  in  terms  of 
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(t^;,  x,..,  x;|!,  x.;!)  by  equations  similar  to  (IV-16),  the  metric  (IV-21)  becomes 
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(IV -23) 


ds 


2 


2,,  2 


c  dt 


r2  (de|  +  sin2  Q^dtp^)  . 


The  metric  given  by  (IV-21)  or  (1V-23)  is  called  the  Schwarzschild  exterior  solution  of  the  Einstein 
field  equations.  The  solution  is  not  valid  inside  the  body,  so  that  the  apparent  singularity  in  the 
metric  for  r^  =  2a  does  not  really  exist. 

If  we  write  out  (IV- 14)  for  the  motion  of  a  body  of  small  mass  in  the  metric  (IV-21),  we  will 

Jr 

obtain  Newton's  equations  of  motion  with  a  small  correction  R  on  the  right-hand  side.  We  then 
suppose  that  this  same  correction  applies  to  the  Newtonian  equations  of  motion  of  a  planet  with 
non-negligible  mass  acted  on  by  the  gravitational  attraction  of  the  Sun  and  other  planets  and  by 
other  forces,  obtaining  tIII-4)..  Of  course,  the  rigorously  correct  procedure  would  be  to  derive 
the  equations  of  motion  in  a  completely  relativistic  manner,  with  the  equations  for  the  comparison 
of  theory  and  observation  also  being  derived  according  to  the  general  theory  of  relativity.  How¬ 
ever,  given  the  limitations  stated  in  Sec.  I,  we  continue  with  the  less  rigorous  procedure  of  using 
the  relativistic  equations  of  motion  of  a  planet  of  small  mass  in  the  gravitational  field  of  the  Sun 

to  correct  the  Newtonian  equations  of  motion  of  a  planet. 
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The  coordinate  system  (t^,  x^,  x#,  x^)  could  be  changed  very  slightly  and  the  equations  of 

motion  of  a  particle  of  small  mass  would  still  have  the  appearance  of  the  Newtonian  equations  of 

motion  with  a  small  (but  different)  correction  on  the  right-hand  side.  If  we  are  to  follow  the  plan 

of  correcting  the  Newtonian  equations  of  motion,  this  correction  should  be  obtained  in  the  general 

relativistic  coordinate  system  which  most  closely  fits  the  Newtonian  coordinate  system.  The 
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only  reason  that  the  coordinate  system  (tJic,  x*,  x*,  x,,)  with  metric  (IV-21)  could  be  this  "best- 

fitting"  coordinate  system  is  the  apparent  simplicity  of  the  metric  (IV-21). 
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In  a  coordinate  system  (t,  x  ,x  ,x  )  with  flat  metric  (IV-15),  the  d'Alemfcertian  operator  □ 
is  defined  on  a  function  f  by  the  equation 


?  2  2  2 

1  9f  9^f  9  f  9  f 

2  „,2  1,2  2,2  3,2 

c  9t  9(x  )  9(x  )  9(x  ) 


(IV-24) 


Note  that  Dt  =  0,  Dx^  =  0  (j  =  1,  2,  3).  The  natural  differential  geometric  generalization  of  the 
d'Alembertian  operator  to  a  coordinate  system  (x°,  x1,  x2,  x3)  with  metric  ds2  =  g  dx*1  ®  dx^  is 
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(IV-25) 
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where  the  summation  on  the  Greek  indices  runs  over  0,  1,  2,  3.'  A  coordinate  system  (x  ,  x  ,  x  ,  x  ) 

is  harmonic  if  Ox*1  =  0  (p  =  0,  1,  2,  3).  Given  any  coordinate  system  (x£,  x*,  x  2,  x3)  and  metric, 
a  new  coordinate  system  (x^,  x*,  x2,  x3)  which  is  harmonic  can,  in  general,  be  found.  (We 
would  have  to  find  four  independent  solutions  of  the  linear  hyperbolic  partial  differential  equa¬ 
tion  with  nonconstant  coefficients  Of  =  0,  which  can  be  done  because  Cauchy's  problem  can 
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be  solved  for  this  type  of  partial  differential  equation.  )  Suppose  we  are  concerned  with  an 

insular  distribution  of  matter  contained  in  a  coordinate  neighborhood.  We  may  assume  that 

this  coordinate  neighborhood  extends  off  to  infinity,  with  the  insular  system  being  contained  in 

28 

a  spatially  bounded  part  of  the  coordinate  neighborhood.  Then,  Fock  proves  that  if  certain 
natural  conditions  are  satisfied  by  the  metric  at  spatial  infinity,  a  harmonic  coordinate  system 
on  the  coordinate  neighborhood  is  defined  uniquely  up  to  a  Lorentz  transformation.  These  con¬ 
ditions  essentially  state  that  the  metric  goes  sufficiently  fast  to  the  flat  space  value  (III- 15)  at 


tSee,  for  example,  Ref. 26. 
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spatial  infinity,  and  that  no  gravitational  waves  impinge  on  the  insular  system  from  the  outside. 

In  the  case  of  the  isolated  spherically  symmetric  body  considered  at  the  beginning  of  this 

12  3 

section,  the  harmonic  coordinate  system  (t,  x  ,  x  ,  x  )  can  be  made  unique  up  to  a  rotation  of  the 
spatial  axes  by  specifying  that  the  world  line  followed  by  the  center  of  the  spherically  symmetric 

body  is  given  by  x1  =  0  (i  =  1,  2,  3).  In  these  harmonic  coordinates,  the  metiic  tensor  has  com- 

.  29 
ponents 

,r  —  a  k  2  A 

^oo  r  +  a  c  ’  ®oi  " 


g..  =  -(l  +  “)z  6..-(L±°L)  {£)2 
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(IV-26) 


where  a  is  given  by  (IV-22).  In  the  spherical  coordinate  system  (t,  r,  0,  tp),  defined  in  terms 
of  (t,  x1,  x2,  x3)  by  (IV-16),  we  have2^ 


ds2  =  (jhjrf)  c2dt2-  dr2  -  Ha)2  (d©2  +  sin2  ©d?2)  .  (IV-27) 

The  relation  between  the  (t#,  r*.,  0^,  <p coordinate  system  of  (IV-23)  and  the(t,  r,  0,  <p)  coordi¬ 
nate  system  of  (IV -27)  is  obviously  given  by 


t...  =  t  ,  ry.  =  r  +  a 

v  v 


0*  =  0  -  =  <P  (IV-28) 
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so  that  the  transformation  between  the  (t#,  x,,.,  xif,  x*)  coordinate  system  of  (IV-21)  and  the 
12  3 

(t,  x  ,  x  ,  x  )  coordinate  system  of  (IV-26)  is  given  by 


■t*  =  t  ,  x*  =  x1  (1  +  2L)  ,  i  =  1,  2,  3 


(IV-29) 


The  Schwarzschild  metric  has  also  been  expressed  in  what  are  called  isotropic  coordinates. 

-  _1  2  _3 

In  isotropic  rectangular  coordinates  (t,  x  ,  x  ,x  ),  the  Schwai  zschild  metric  has  components 
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while  in  isotropic  spherical  coordinates  (t,  r,  0,  <p),  it  is  given  by 

ds2  =  — - c2dt2  -  (1  +  =^)4  (dr2  +  r2d©2  +  r2  sin  2  ©d?2) 

(1  +  #> 

(see  Ref,  30).  Comparing  (IV-23),  (IV-28)  and  (IV-31),  we  see  that 


(IV-30) 


(IV-31) 
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which  implies  that 
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(IV— 32) 


(IV-33) 


There  are,  of  course,  infinitely  many  other  coordinate  systems  besides  these  three  in  which 
the  Schwarzschild  metric  can  be  expressed,  but  these  three  are  generally  used  in  the  literature. 
We  shall  now  derive  the  tquations  of  motion  of  a  particle  of  small  mass  in  each  of  the  coordinate 
systems  discussed  above,  even  though  we  have  reason  to  believe  that  the  harmonic  coordinates 
are  closest  to  the  Newtonian  coordinates. 

First,  we  note  that,  by  Ref.  29, 
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by  Ref.  24, 
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and,  by  (IV-30), 
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Thus,  by  definition  (IV-4),  we  have  in  the  (t,  x  ,x  ,x  )  coordinate  system  that 
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(IV-37) 


3xJ  3x  3x 


wiih  exactly  similar  equations  holding  in  the  starred  and  barred  coordinate  systems.  Equation 
(IV-14)  for  a  geodesic  then  becomes 

d  t  _  y  p  o  dt  dx^  _  q 
ds2  oj  ds  ds  " 


Since 
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the  second  equation  in  (IV-38)  can  be  written 


(1V-38) 
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Now,  by  the  first  equation  in  (IV-38),  we  have 
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so  that  we  can  finally  write 
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(IV-39) 


Exactly  similar  equations  are  valid  in  the  starred  and  barred  coordinate  systems. 
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Using  (IV-26),  (IV-34)  and  (IV-37),  v.re  perform  a  simple  calculation  in  the  (t,  x  ,x  ,  x  ) 
coordinate  system  that  gives 
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(IV -40) 


Similarly,  using  (IV-21),  (IV-35)  and  (IV-37),  we  find  that  in  the  (t#,  x.*,  x2,  x3)  coordinate 
system 
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Finally,  using  (IV-30),  (IV-36)  and  (IV-37),  we  find  that  in  the  (t,  x  ,  x  ,x  )  coordinate  system 
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(IV -42) 


Formulas  (IV-39)  and  (IV-40)  show  that  the  equations  of  motion  of  a  small  mass  m  the 
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(t,  x  ,  x  ,  x  )  coordinate  system  are 
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Similarly,  formulas  (IV-39)  and  (IV-41)  show  that  the  equations  of  motion  of  a  particle  of  small 
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mass  in  the  (t^,  x*,  x^,  x^)  coordinate  system  are 
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Finally,  formulas  (IV-39)  and  (IV-42)  show  that  the  equations  of  motion  of  a  particle  of  small 
-  _1  _2  _3 

mass  in  the  (t,  x  ,  x  ,  x  )  coordinate  system  a  re 
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In  the  case  of  Mercury, 
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Thus,  using  the  fact  that  for  small  z 
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we  drop  all  terms  in  (IV-44),  (IV-46)  and  (IV-48)  which  contain  (a/r)  or  (a/r)  (v/c)  as  factors, 
obtaining. 
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(IV-51) 


Equations  (IV-43),  (IV-45)  and  (IV-47)  are  invariant  under  rotation  of  the  coordinate  axes. 
From  this  it  easily  follows  that  the  motion  given  by  these  equations  lies  in  a  plane  in  each  of  the 
three  coordinate  systems. 

By  coordinate  transformations  (IV-28),  (IV-29),  (IV-32)  and  (IV-33),  the  curves  in  the  three 
coordinate  systems  given  by  {IV-43),  (IV-45)  and  (IV-47)  are  exactly  similar,  even  though  these 
equations  have  dissimilar  appearance.  As  is  well  known,  these  curves  are  ellipses  with  advanc- 

o  i  on  -Jo 

ing  perihelions.  '  '  "  The  expressions  for  the  periods  of  these  ellipses,  in  terms  of  the  semi¬ 

major  axes  of  the  ellipses,  will  vary  in  the  different  coordinate  systems  because  of  relations 
(IV-28)  and  (IV-32). 

We  have  two  candidates  for  the  relativity  term  in  (III-4)  and  (III-12).  Because  of  the  har¬ 
monic  and  isotropic  criteria,  and  for  the  sake  of  definitiveness,  we  choose  (IV-49)  [or  equiva¬ 
lently  (IV-51)]  to  be  this  term.  In  the  notation  of  Sec.III-A.,  it  is 
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(TV-52) 


where  we  have  multiplied  (IV-49)  by  a  dimensionless  constant  R^..  If  we  perform  a  least-squares 
analysis  on  the  value  of  Rf  and  other  parameters  to  fit  theory  and  observation,  Rf  will  converge 
to  the  value  1  if  the  relativity  correction  belongs  in  the  equations  of  motion,  or  to  the  value  0 
if  the  Newtonian  theory  is  correct.  By  (IV-22), 

yM 

a  =  — Y  (IV-53) 

c 


where  c  is  the  velocity  of  light  at  a  large  spatial  distance  from  the  Sun. 


C.  METHOD  OF  SOLVING  THE  PROBLEM  OF  THE  MOTION  OF  A  SYSTEM 
OF  MASSES  IN  GENERAL  RELATIVITY 


If  we  raise  indices  in  the  Einstein  field  equations  (IV-7),  we  obtain 

RH"  _  i  g**"  R  =  -  kT^  .  (IV-F  !) 

These  equations  are  nonlinear  and  hyperbolic  in  the  unknown  functions  g^ .  The  fact  that  they 
are  hyperbolic  implies  that  gravitational  waves  can  exist.  Their  nonlinearity  allows  us  not  only 
to  determine  the  potential  g^u,  but  also  the  mass  tensor  T^P,  i.e.,  the  motion  of  the  masses. 

In  all  other  field  theories,  such  as  Maxwell's  for  the  electromagnetic  field  or  Newton's  for  the 
gravitational  field,  the  field  equations  are  linear  and  the  equations  for  the  motion  of  bodies  in 
the  field  are  separate  from  and  additional  to  the  field  equations.  But  in  Emstein's  theory,  the 
equations  of  motion  are  contained  in  the  equations  for  the  field. 

The  divergence  of  the  left  side  of  the  Einstein  field  equations  vanishes. 
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(IV -55) 
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(This  was  one  of  the  attributes  which  led  Einstein  to  choose  this  tensor  for  the  left  side  of  his 
equation.)  Thus,  we  have 


v  +  T<7Vr  11  +  T^r  u  =  0  .  (IV-56) 

^  3x^  ^ 

The  simultaneous  solution  of  (IV -54)  and  (IV-56)  will  determine  the  field  and  the  motion  of  the 
masses;  of  course,  (IV-56)  is  a  consequence  of  (IV-54). 

0  "1  2  3 

To  obtain  an  approximate  expression  for  the  mass  tensor,  let  p(x  ,  x  ,  x  ,  xJ)  be  (in  some 
sense)  the  invariant  density  of  matter  in  the  space -time  universe.  We  suppose  that 


= 


(IV-57) 


where  s  is  the  proper  time  of  the  element  of  matter  at  the  point  (x^,  x\  x^,  x^).  If  the  element 
of  matter  is  following  a  world  line  x^1  =  x^(s),  then  the  defining  property  of  s  is 
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If  we  imagine  that  we  are  concerned  with  a  particle  of  small  mass  and  small  dimensions  which 

has  negligible  effect  on  the  gravitational  field,  (IV-56)  and  (IV-57)  imply  that  the  particle  follows 
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a  time-like  geodesic  through  the  space-time  universe.  Thus,  assumption  (IV-14)  is  not  really 

an  assumption,  but  is  a  consequence  of  the  Einstein  field  equations. 

The  method  of  solving  the  field  equations  for  the  field  and  the  motion  of  the  masses  is  pre- 
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sented  by  Fock,  and  by  Infeld  and  Plebanski.  The  latter  follow  the  work  of  Einstein,  Infeld 
36 

and  Hoffmann  and  assume  that  the  masses  are  point  singularities  of  the  field,  so  that  the  mass 
tensor  is  zero  everywhere  except  along  the  world  lines  of  the  particles,  where  it  is  given  by 
delta  functions.  Fock  assumes  a  continuous  distribution  of  matter  concentrated  in  a  finite  num¬ 
ber  of  regions,  so  that  the  mass  tensor  is  differentiable  everywhere,  and  is  zero  outside  of  the 
finite  number  of  regions.  The  methods  used  by  Fock  and  by  Infeld  and  Plebanski  are  approxi¬ 
mation  procedures  and  are  essentially  the  same;  Fock  assumes  that  he  is  always  working  in  a 
harmonic  coordinate  system,  while  Infeld  and  Plebanski  make  supplementary  coordinate  condi¬ 
tions  at  each  step  in  the  approximation. 

To  be  specific,  let  us  outline  Fock's  procedure  with  a  continuous  distribution  of  matter 
concentrated  in  a  finite  number  of  regions.  We  first  assume  expression  (IV-57)  for  the  mass 
tensor  T (at  a  later  stage  in  the  approximation,  we  can  assume  a  more  sophisticated  form  of 
the  mass  tensor  using  the  fact  that  the  bodies  are  elastic).  Then  we  solve  (IV-54)  for  the  gravi¬ 
tational  potential  g1^  to  first  order  in  v/c,  obtaining  (IV-19)  with  some  additional  terms  of  the 
form  dt  ®  dx1  times  quantities  involving  the  velocity  of  the  matter  generating  the  field.  This  is 
called  the  Newtonian  approximation.  We  use  this  solution  for  the  potential  to  write  equations 
(IV-56)  for  the  mass  tensor  T1^  to  first  order  in  v/c.  The  solution  of  these  equations  is  used 
to  solve  (IV-54)  for  the  gravitational  potential  g1^  to  second  order  in  v/c.  This  solution  for  the 
potential  is  then  used  to  write  equations  (IV-56)  for  the  mass  tensor  T^  to  second  order  in  v/c. 
We  could,  in  principle,  continue  this  procedure  indefinitely,  but  the  solutions  and  equations  in 
this  post-Newtonian  approximation  are  accurate  enough  for  our  purposes. 
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At  whatever  stage  we  stop  in  the  approximation  procedure,  we  will  have  found  in  a  specific 

coordinate  system  a  system  of  second  order  ordinary  differential  equations  for  the  motions  of 

the  masses,  and  an  expression  for  the  gravitational  potential  in  terms  of  the  motions  of  masses. 

12  3 

If  the  coordinate  system  is  it,  x  ,  x  ,  x  ),  t  can  be  made  the  independent  variable  of  the  equations 
of  motion.  Numerical  integration  of  the  equations  of  motion  will  determine  ephemerides  in  this 
specific  coordinate  system  of  the  planets  as  functions  of  t.  The  relation  between  the  coordinate 
time  t  and  the  proper  time  r  of  an  atomic  clock  on  the  surface  of  the  Earth  following  a  world 
line  y(s)  (s  =  t)  is  given  by  (IV-12).  This  integral  can  be  evaluated  because  we  know  the  gravi¬ 
tational  potential  in  our  specific  coordinate  system.  Knowing  the  general  relativistic  theory  of 
radar  and  optical  observations  of  a  planet,  we  can  comp,  te  the  theoretical  values  of  observations 
made  at  given  instants  of  atomic  time.  Then,  making  a  least-squares  adjustment  to  the  initial 
conditions  and  parameters  appearing  in  the  theory  of  motion,  we  can  determine  the  general 
relativistic  ephemerides  which  best  fit  observation. 

This  report  is  concerned  with  Newtonian  theory  and  any  general  relativistic  corrections 

that  are  easily  obtained.  The  procedure  outlined  above  can  and  should  be  documented  and 

21  35 

pursued;  this  we  hope  to  do,  following  Fock  and  Infeld  and  Plebanski. 
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V.  EQUATIONS  FOR  PARTIAL  DERIVATIVES  OF  POSITION  AND  VELOCITY 
WITH  RESPECT  TO  INITIAL  CONDITIONS  AND  PARAMETERS 


A.  PLANET  CASE 


Let  (x 


ps' 


. ,  x  )  denote  the  components  of  position  and  velocity  of  a  planet  relative  to  the 
ps 


Sun.  Equation  (III-4)  for  the  motion  of  the  planet  can  be  written  in  the  form 
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(V-l) 


Let  a  be  a  parameter  upon  which  the  motion  of  the  planet  depends,  such  as  an  initial  condition, 
a  planetary  mass,  the  second  harmonic  of  tne  Sun,  etc.  Differentiating  system  (V-l)  with  respect 
to  a,  we  see  that  the  quantities  (3x  1  /3a,  . . . ,  9x  ^  /3a)  satisfy  the  differential  equations  system 
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(V-2) 
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We  have  3xZ0/9a  =  0  (t  =  1, .  . . ,  6)  for  all  parameters  a  which  are  not  initial  conditions. 

Opb 

The  value  of  /9a  for  an  initial  condition  a  depends  on  the  specific  initial  condition.  For 

-j  ^ 

example,  if  we  let  a  =  x"pg,  then  9x0pS/3«  =  <5^.  We  shall  choose  the  initial  conditions  as 

(0p,  •  •  ■  .  )  =  (a,  e,  i,  n,  w,fo),  the  orbital  elements  of  the  elliptic  orbit  osculating  to  the  true 

orbit  of  the  planet  at  the  initial  time  t  ,  because  we  will  be  integrating  differential  equations 
k  /  i  ^ 

system  (V-2)  for  3x  /3/?J  (j,  k  =  1,  .  .  . ,  6)  by  means  of  Encke's  method  as  explained  in  Sec.  VI-A, 
ps  p 

Also,  we  are  going  to  use  the  results  of  the  integration  to  make  a  least -squares  correction  to 

the  initial  conditions,  and  it  is  more  meaningful  physically  to  adjust  (/?*,  .  .  ,  than  to  adjust 

16  1  6  1  6"  P 

(xQps,  .  .  .  ,xQps).  The  relation  between  (/3p ,  . .  . ,  /?p )  and  (xQps,  . .  •  ,xopp/  is  given  in  Secs.  II -B 

and  II-C,  while  the  values  of  3xk  /3 (3^  are  given  in  formulas  (11-28)  to  (11-33)  with  t  =  t  . 

ops  p  o 
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If  we  desire  to  take  account  of  a  possible  time  variation  of  the  gravitational  constant,  we 
might  suppose  that 


YMs  =  (yMs)o  [1  +  X(t-to)] 


( V-3) 


where  X  is  a  parameter  to  be  determined  by  comparing  theory  with  observation.  We  then  have 


9(YMS) 


=  (TMS)0  (t-tc) 


(V-4) 


with  the  3(yMs)/3a  term  in  (V-2)  being  zero  for  all  other  parameters  a.t 

The  term  3{M  /M  )/da  in  (V-2)  is  zero  for  all  parameters  a,  except  for  «_•  =  M  /M  ,  in 
p  s  p  s 

which  case  it  takes  the  value  1. 

To  determine  the  term  3 ft  /3a  in  (V-2),  we  differentiate  (IJI-3)  with  respect  to  a,  obtaining 
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(V-5) 


if  a  is  not  the  mass  of  a  perturbing  planet.  Here,  we  have  to  assume  that 


3xk 

-Jl  -  0 

3a 


3x.k  3x.k 

JP  =  ]S 

3a  da 


3xk 
PS  _ 
3a 
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(V-6) 
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since  it  is  supposed  that  (x.  ,  x.  ,  x.  )  are  given  as  definite  functions  of  time.  If  a  is  the  mass 

Is  3s 

of  some  perturbing  planet,  a  =  M,/M0,  then  we  must  add  the  following  term  to  (V-5): 

1  s 


/xk  xk 

\r.  r. 
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k  =  1,  2,  3 


(V-7) 


If  a  =  X,  we  must  add  this  term  to  (V-5): 


■  a  s  \r.  r. 
3=1  '  JP  js 


k  =1,  2,  3 


(V-8) 


If  a  is  an  initial  condition  for  planet  p,  assum;tion  (V-6)  that  3x.  /3a  is  zero  relative  to 

*  JS 

3xDg/3a  is  certainly  justified.  However,  if  a  is  some  parameter  other  than  an  initial  condition, 

3x.k/3a  could  be  comparable  or  even  larger  than  3x  Jda  so  that  expression  (V-5)  for  3ft  /3a 
js  ps 

would  be  incorrect  in  nature.  For  use  in  PEP,  (V-5)  is  exactly  correct  when  the  perturbing 
planet  input  magnetic  tape  used  in  evaluating  (III-3)  is  kept  the  same  between  iterations,  but  not 
when  the  ephemerides  on  the  perturbing  planet  tape  are  replaced  by  the  results  of  just  completed 


t  We  do  not  consider  (yMs)0  as  a  parameter  to  be  adjusted  because  it  is  the  usual  practice  in  celestial  mechanics 
to  set  i/(yMs)0  =  0.01720209895,  which  defines  the  unit  of  length  (the  Astronomical  Unit)  once  the  unit  of  time 
has  been  specified. 
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integrations  for  use  in  the  next  iteration.  In  any  case,  because  (V-5)  contains  the  factor  {M./M  ), 
it  is  less  important  than  (V-7)  or  the  first  term  on  the  right  of  the  second  equation  in  (V-2) J  S 
Since  we  only  need  approximate  values  for  the  partial  derivatives  in  the  iterative  process  of 
finding  least-squares  corrections  to  the  various  parameters,  assumption  (V-6)  is  thus  reasonable. 
From  an  operational  standpoint,  it  is  necessary. 

It  is  probably  sufficiently  accurate  to  suppose  that  9Rk/9a  =  0  and  9Sk/9a  =  0,  except  v/hen 
a  =  Rf  or  a  =  S,/Rs,  respectively.  However,  for  the  sake  of  completeness,  we  shall  evaluate 
these  quantities.  First,  differentiating  (IV-52)  with  respect  to  5,  we  see  that 
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We  have  denoted  the  parameter  with  which  we  differentiate  (IV-52)  by  a  so  that  there  will  be  no 
confusion  with  the  gravitational  radius  of  the  Sun  a  appearing  in  this  formula.  If  a  =  Rj,  we 
must  add  the  following  term  to  (V-9): 
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If  a  =  X,  we  must  add  this  term  to  (V-9) 
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(V-12) 


Second,  differentiating  (III-50)  with  respect  to  a,  we  see  that 
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If  a  =  S^/R  ,  we  must  add  the  following  term  to  (V-13): 
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If  a  =  Mp/Mg,  we  must  add  this  term  to  (V-13): 


yMs  / Rs \ 2  S2  [xds  ,15  2_3 
r 2  kps'  R2  [^ps  2  2>  g  3k  ' 


k  =  1,  2,  3 


If  a  =  X,  we  must  add  this  term  to  (V-13): 
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B.  EARTH-MOON  CASE 

Let  *xcs’  •  •  •  ’  XCS>  and  *xme'  * '  •  *  xme^  denote  the  components  of  position  and  velocity  of  the 
Earth-Moon  barycenter  relative  to  the  Sun,  and  of  the  Moon  relative  to  the  Earth,  respectively. 
These  components  will  be  considered  as  primary  quantities  determined  by  integrating  the  equa¬ 
tions  of  motion.  The  components  of  position  and  velocity  of  the  Earth  and  Moon  relative  to  the 

Sun  ^xes’  •  ’  * '  xes^  and  *xms'  ■  •  • '  xms^  are  determined  from  them  by  formulas  (III-6).  Equations 
(III -12)  for  the  motions  of  the  Earth-Moon  barycenter  and  the  Moon  can  be  written  in  the  form 
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xcks  =  xokcs  '  xcks3  =  xokcs  when  t  =  tQ 
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Let  a  be  some  parameter  upon  which  the  motions  of  the  Earth-Moon  barycenter  and  of  the 

Moon  depend.  Differentiating  systems  (V-18)  and  (V-19;  with  respect  to  a,  we  see  that  the 
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Let  (0  ,  . . .,  0  )  and  (0  ,  ....  0  )  be  the  orbital  elements  of  the  elliptic  oi  bits  osculating 

at  the  initial  time  to  the  true  orbits  of  the  Earth-Moon  barycenter  around  the  Sun  and  of  the 

Moon  around  the  Earth,  respectively.  Since  we  shall  be  integrating  the  differential  equations 

for  the  partial  derivatives  with  respect  to  initial  conditions  by  means  of  Encke's  method,  we 

choose  the  initial  conditions  with  respect  to  which  we  take  partial  derivatives  to  be  (0^,  . . . ,  p^) 

and  (0  * ,  . . . ,  P t).  The  initial  conditions  9x^/9 p\  and  dx^/dp^  (j,  k  =  1,  . . . ,  6)  are  then 
m  m  ocs  c  oiuc  m 

determined  by  the  elliptic  orbit  formulas  of  Sec.II-D.  Of  course,  we  have 


j,  k  =  1,  ....  6 


( V-22) 


Further,  the  initial  conditions  of  (V-20)  and  (V-21)  are  zero  for  any  parameter  a  not  an  initial 
condition. 

We  divide  the  initial  conditions  and  parameters  appearing  in  the  theories  of  motion  of  the 
Earth-Moon  barycenter  and  of  the  Moon  into  the  following  three  classes: 


Initial  osculating  orbital  elements  of  orbit 
of  Earth-Moon  barycenter  about  the  Sun 
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Mass  of  j  perturbing  planet 
Mass  of  Sun 
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Second  harmonic  of  the  Sun 
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Initial  osculating  orbital  elements  of  orbit 
of  the  Moon  about  the  Earth 
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Mass  of  Earth  +  Moon 
Mass  of  Sun 
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Mass  of  Moon 
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(V-25) 


Time  variation  factor  for  gravitational  ^  | 

constant  J 

For  parameters  a  of  the  form  (V-23),  we  shall  assume  that  9x^/9  a  =  0,  and  for  (V-24)  that 
.  me 

9xK  /da  =  0.  We  make  no  such  assumptions  concerning  parameters  (V-25).  Comparing  the 
CS 

magnitude  of  the  A  term  in  (III- 13)  given  by  (III-32)  with  the  magnitude  of  the  perturbing  planet 

effects  in  Table  II,  we  see  that  it  is  indeed  reasonable  to  assume  that  9x^/9 pj^  =  0,  since  we 

must  assume  that  <f>  is  independent  of  the  initial  conditions  for  the  perturbing  planets.  The 

assumption  that  9x^  /90^  =  0  is  not  precisely  true,  but  according  to  Table  III  these  derivatives 
me  c  k  /  i 

are  much  smaller  than  the  anc*  we  have  to  make  some  such  assumption  to  make  our 

problem  manageable.  Comparing  Tables  II  and  III,  we  see  that  it  is  reasonable  to  assume  that 
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the  3xk  /3M.  are  zero  relative  to  the  3xk  /3M.  .  (We  repeat  that  in  the  iterative  determination 
me'  js  cs'  js  1 

of  the  least-squares  correction  to  the  various  parameters,  only  the  approximately  exact  values 
of  the  partial  derivatives  with  respect  to  these  parameters  are  needed,  not  the  exact  values.) 

For  parameters  of  the  type  (V-23)  for  which  we  assume  9xr  ie/3o:  =  0,  the  equations  in  (III-6) 


give 


Sx" 

es 

3a 


9x 


cs 


3x 


3a 

k 


ms 


3x 


cs 


3a 


3a 


k  =  1,2,3 


(V-26) 


We  now  determine  the  quantities  on  the  right  side  of  (V-20)  for  these  parameters  a.  First,  of 
course,  3(yMs)/3a  =  0  and  (3/ 3a)  (Mc/Mg)  =  0.  Next,  we  have 
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The  expressions  for  3Rk/3a  and  3Sk/3a  are  given  in  formulas  (V-9)  through  (V-17),  with  p  re¬ 
placed  by  c.  Finally,  differentiating  (HI-9)  with  respect  to  a,  we  see  that 
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where  we  have  made  assumptions  analogous  to  those  of  (V-6),  namely, 
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If  a  =  M./M  ,  we  must  add  the  following  term  to  expression  (V-28)  for  3$k/3a: 
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For  the  initial  conditions  (0*  . . . ,  4)  of  (V-24),  we  assume  that  3x^/3/?  ^  =  0.  By  (III-6), 
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We  now  determine  the  quantities  on  the  right  of  (V-21)  for  a  =  4’  First,  of  course, 

3(yM  )/3/4  =  0  and  { B/ 3/3  ^  )  (M  /M  )  =  0.  Next,  differentiating  (III- 11)  with  respect  to  /4  and 
s  m  m  c  s  m 

using  (V-31),  we  see  that 

■^  =  VM  fs  y  — 52  ('5a  ^ 


3B  „  V  ine  r  m  es‘ es 

— r =  yMs  3  2i  — r  TvT"  •  t 

3/3'*  S  ,  ,  3/3 J  \Mc  r 
Fm  *■  i  =  l  Fm  '  es 

3xk  /M  .  IVI  ,  \1 

_  me  I _ m  1  ,  _ e  1  \ 

„„j  \M  3  M  3  / 

30 J  \  c  r  „  cr  /  , 

•  '  ac  me*  •* 


k  =  1,  2,  3 


(V-32) 


5-^2  C  3  .2 


9xme 

30^ 

Fm 


IVI  .  M.  , 

m  1  ,  e  l 

M  3  +  M  3 

c  r.  c  r. 
le  1m 


k  t  k  i 

M  x.  x.  M  x.  x. 

m  le  le  e  im  1m 

M  5  M  5 

c  r.  c  r. 

le  im 


k  =  1,  2,  3 


In  (V-33)  we  have  assumed  that 


9xis 

=  0 

< 

3x.k 

le 

ax.1! 

IS 

axk 

es 

M 

m 

q  k 

3x 

me 

< 

30;* 

m 

30;* 

m 

Mc 

9^ 

0xim 

q  k 
3x. 
is 

q  k 

9xms  _ 

M 

9xk 

e  _ me 

30-* 

m 

3/3  ■* 
m 

s4 " 

"M 

c  30;* 
m 

k  =  1,  2,  3 


(V-33) 


(V-34) 


Finally,  differentiating  (III-79)  with  respect  to  0^,  we  see  that 
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where,  by  (III-80), 
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We  now  consider  the  parameters  M  =  M  /M  and  M  =  M  /M  .  Regarding  these  as 
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independent  parameters,  formulas  (III-6)  imply  that 
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The  quantities  on  the  right-hand  sides  of  (V-20)  and  (V-21)  for  a  =  Mcg  are  as  follows.  First, 


of  course,  a(yMg)/8Mcg  =  0  and  (8/dMcg)  (Mc/M  )  =  1.  Next,  we  have 
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The  expressions  for  3Rk/aMcs  and  3Sk/aMcs  are  given  in  formulas  (V-9)  through  (V-17),  with 
p  replaced  by  c.  Differentiating  (111-9)  with  respect  to  Mcg,  we  see  that 
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where  we  have  assumed  that 
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Differentiating  (III-ll)  with  respect  to  M  ,  we  see  that 
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where  we  have  used  (V-44)  in  deriving  (V-46).  Let  denote  the  right  side  of  (V-35),  with 

partial  derivatives  with  respect  to  /3^  replaced  by  partial  derivatives  with  respect  to  Mcp. 

Then,  differentiating  (III-79)  with  respect  to  M„_,  we  see  that 
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The  quantities  on  the  right-hand  sides  of  (V-20)  and  (V-21)  for  a  =  M  ,  are  as  follows. 
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The  expressions  for  9Rk/9Mmc  and  9Sk/9Mmc  are  given  in  formulas  (V-9)  through  (V-17),  with 
p  replaced  by  c.  Differentiating  (III-9)  with  respect  to  Mmc,  we  see  that 
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where  we  have  assumed  that 
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Differentiating  (III-ll)  with  respect  to  Mmc,  we  see  that 
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where  we  have  used  (V-50)  in  deriving  (V-52).  Differentiating  (III-79)  with  respect  to  M  ,  u 
see  that  3H /3M  c  is  given  by  (V-35),  with  partial  derivatives  with  respect  to  0^  being  re¬ 
placed  by  partial  derivatives  with  respect  to  Mm  . 

We  now  consider  the  time  variation  parameter  X  in  the  gravitational  constant  as  given  in 
(V-3).  Formulas  (III-6)  imply  that 
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The  quantities  on  the  right  sides  of  (V-20)  and  (V-21)  for  a  -  X  are  as  follows.  First,  of  course, 

3(yM  )/ 3X  is  given  by(V-4)  and(3/3X)(M  /M  )  =  0.  Next,  we  have  that(9/3X)  [(M  /M  )  (xk  / r3  ) 

(M  /M  )  (xK  /r  )]  is  given  by  the  right  side  of  (V-42)  with  partial  derivatives  with  respect  to 
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M  being  replaced  by  partial  derivatives  with  respect  to  X.  The  expressions  for  3RK/3x  and 
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If  we  let  B^  and  \I>^  denote  the  right-hand  sides  of  (V-45)  and  (V-46),  respectively,  with  the 
partial  derivatives  with  respect  to  M£S  being  replaced  by  partial  derivatives  with  respect  to  X, 
differentiation  of  (III-ll)  gives 
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If  we  let  H.  denote  the  right-hand  side  of  (V-35),  with  partial  derivatives  with  respect  to  /3 J 
A  m 

being  replaced  by  partial  derivatives  with  respect  to  X,  differentiation  of  (III-79)  gives 
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We  have  not  derived  the  differential  equations  satisfied  by  the  partial  derivatives  of  the  posi¬ 
tion  and  velocity  of  the  Moon  with  respect  to  the  higher  harmonics  of  the  gravitational  potentials 
of  the  Earth  and  Moon,  because  the  gravitational  potential  of  the  Earth  has  been  determined  quite 
accurately  from  the  motion  of  artificial  Earth  satellites,  and  the  gravitational  potential  of  the 
Moon  will  be  determined  quite  accurately  in  the  near  future  from  the  motion  of  artificial  lunar 
satellites. 
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VI.  ENCKE'S  EQUATIONS  OF  MOTION  AND  ENCKE’S  EQUATIONS 
FOR  PARTIAL  DERIVATIVES  WITH  RESPECT  TO  INITIAL  CONDITIONS 


A.  PLANET  CASE 

Let  (x  1  ,  . .  . ,  x  ^  )  denote  the  components  of  position  and  velocity  of  a  planet  relative  to  the 
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Sun,  satisfying  (V-l)  with  initial  conditions  (xopg,  ■  •  -'x0pS)  anc®  with  the  gravitational  constant 

yM_  being  given  by  (V-3).  Let  (y  *  . . . ,  y  ^  )  be  the  solutions  of  the  equations  in  (11-63)  with 
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p  :=  (yMs)o  (1  +  (Mp/Mg)],  and  with  initial  conditions  (xopg,  •  •  . ,  xQps).  The  quantities 

(y  i  . . . ,  y  ^  )  are  the  components  of  position  and  velocity  in  the  elliptic  orbit  osculavmg  to  the 
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i  +  _J E 
l  M 

'  S' 


\  [/  1  1  \  k  ?ps 

Wns  ns'  pS 


k  =  1,  2,  3 


(VI-2) 


-X(t  -  tQ) 


■f  fik  +  Rk  +  Sk 


when  t  =  t 


I  1  2  2  2  3  2  1  6 

where  ppg  =  J(yps)  +  (yps)  +  (yps)  .  The  quantities  (ypg,  . .  . ,  ypg)  are  known  as  functions 

of  time  from  the  formulas  in  Sec.II-B,  so  that  if  we  numerically  integrate  (VI-2)  to  find 

(|pg,  ....  lpS).  the  position  and  velocity  of  the  planet  (xpg,  . .  . ,  xpg)  can  be  determined  from  (VI-1). 

Let  3x^ ■ /3/3  3  (j,  k  =  1,  .  .  . ,  6)  denote  the  partial  derivatives  of  the  position  and  velocity  of 
Ps  P  1  6 

the  planet  relative  to  the  Sun  with  respect  to  the  initial  osculating  orbital  elements  (/?p ,  .  . . ,  p  )  = 

(a,  e,  i,  S2,  w,  lQ).  These  quantities  satisfy  differential  equations  (V-2)  with  initial  conditions 
3xkpg/3^^  (j,  k  =  1, . . . ,  6).  Let  3yps/3/3^  (j,  k  =  1,  . . . ,  6)  be  the  solutions  of  differential  equa¬ 
tions  (11-64).  with  these  same  initial  conditions  and  with  p  =  (yM  )  [1  +  (M  /M  )).  We  define 

so  p  s 


j,  k  =  1, 


(VI-3) 


Subtracting  (11-64)  from  (V-2),  we  see  that 


■bar** 
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dijk 
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dt 
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PJ 
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*PJ 


=  (rMs)o 


,  M  , 
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„  k  k 
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3 

+  3  £ 
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„  f  .  k  f  k  f  » 
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i?r  It5 
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ps 


ps 

3 

V 
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x<‘  -  V 


'3x 


ps 


3 
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"ps  i=l 


k 

ps 


ax 


ps 
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ps 
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ps  i  =  l 

k  > 
ps 


Wp  ■ 


air 
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9R 

30i  30 


a_  M_ 

30p  ^1V1P 
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?k) 
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k+3  n 

=  77  .  =0 

PJ 


when  t  =  t 


i  jf 

ps'pj 


k  =  1,  2,  3 
j  =  1, ...»  & 


(VI-4) 


The  quantities  9yk  /90^  are  known  as  functions  of  time  from  the  formulas  in  Sec.II-D,  so  that 

Jps'  Hp  k 

if  we  numerically  integrate  (VI-4)  to  find  the  tj  .,  the  partial  derivatives  of  the  position  and 

PJ  k  /  i 

velocity  of  the  planet  with  respect  to  the  initial  osculating  orbital  elements  9XpS/90^  can  be  de¬ 
termined  from  (VI-3). 

k 

If  the  quantities  £  get  too  large  as  time  progresses,  a  new  osculating  elliptic  orbit  can  be 
ps 

chosen  and  the  integration  of  system  (VI-2)  can  commence  again  with  initial  conditions  zero.  If 

we  are  also  integrating  system  (VI-4)  and  desire  to  change  Encke  orbits,  the  following  procedure 

1  6 

must  be  followed.  Let  (0  ,  .  .  .  ,0  )  be  the  osculating  elliptic  orbital  elements  at  the  initial  time 

t  ;  let  (xy*,  .  .  . ,  x* )  and  9xk/90^  (j,  k  =  1,  .  .  . ,  6)  be  the  position,  velocity  and  partial  derivatives 

with  respect  to  initial  conditions  at  the  time  t..c  at  which  we  wish  to  change  Encke  orbits.  These 

quantities  are  known  from  the  numerical  integration  of  the  equations  of  motion  and  the  equations  for 

1  6 

the  partial  derivatives  with  respect  to  initial  conditions  from  time  tQ  to  time  li;. .  Let  (0„. ,  .  .  . ,  0..( ) 
be  the  osculating  elliptic  orbital  elements  at  time  t,,.  determined  from  the  formulas  in  Sec.II-C. 
Integration  of  (VI-2)  and  (VI-4)  with  initial  conditions  zero  from  time  t.,  to  time  t  determines  the 
position  .and  velocity  of  the  planet  (x1,  .  .  .  ,x^)  and  the  partial  derivatives  of  position  and  velocity 
with  respect  to  the  orbital  elements  at  time  ty;,  3xk/90;j  (j,  k  =  i, .  .  . ,  6).  Then,  to  determine  the 
partial  derivatives  with  respect  to  the  orbital  elements  at  time  t  ,  we  must  use  the  relation 


a  k  „  /  „  dx*  90 f' 

*71-  Z[Z  ~i—i 
,  =  1  'i=1 


9x 

90y 


j,k  =  1,  ....  6 


(VI-5) 


D  .  i 

where  the  matrix  30>;./9x).<  is  determined  from  the  formulas  in  Sec.II-E. 

The  elliptic  orbit  position  and  velocity  in  the  new  Encke  orbit  osculating  to  the  true  orbit  at 
time  t#  satisfy  differential  equations  (11-63)  with  p  =  (yMg)^  (1  +  M  /Mg),  where 


(yivy*  =  (yMs)0[l  +X(t*  -to)J 


(VI-6) 


t 
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***•—  —■J*'' 


Thus,  the  factor  (yM  )  in  (VI-2)  with  initial  conditions  zero  at  time  tif  must  be  replaced  by 
s  o 

(yMs)*,  and  the  term  A(t  —  tQ)  replaced  by  a  term  G(t  -  t#)  satisfying 


yM„  =  (yMJ,  (1  +  G(t-t*)] 


(VI-7) 


which  implies  that 


i  +  Mt* -t)  • 


(VI-8) 


Exactly  similar  comments  apply  to  (VI-4). 

B.  EARTH-MOON  CASE 

l  I 

Let  (x  ,  . .  . ,  x  )  and  (x  ,  . . . ,  x  )  denote  the  components  of  position  and  velocity  of  the 
os  os  me  me 

Earth-Moon  barycenter  relative  to  the  Sun  and  of  the  Moon  relative  to  the  Earth,  satisfying 

A  ^  A  ^ 

(V-18)  and  (V-19)  with  initial  conditions  (xocs*  •  •  •  >  xocs)  and  (xome’  •  •  • » xome^  assuirie  that 
the  gravitational  constant  yMg  is  given  by  (V-3).  Let  (ylg,  ... , ,  y^g)  and  (y^,  . . . ,  y^e)  denote 
the  components  of  position  and  velocity  in  the  elliptic  orbits  osculating  at  the  initial  time  to  the 

I  £> 

true  orbits  of  the  Earth-Moon  barycenter  and  of  the  Moon.  The  quantities  (y  ,  ...» ycs)  are  the 
solutions  of  the  equations  in  (11-63)  with  p  =  (yM  )  [1  +  (M  /M  )j  and  with  initial  conditions 
(xocS,  •  ■  • '  X0CS)’  while  the  <luantities  (y^e-  •  •  •  *  y^g)  are  the  solutions  of  (11-63)  with 
H  =  (Yms)0  (mc/ms)  and  with  initial  conditions  (xo*me,  •  ■  • »  xome)-  Let 


=  x  —  y 
cs  J  cs 


c 

:S 

k  = 
k 

me. 


l,  ...»  6 


1 

Subtracting  (11-63)  from  (V-18)  and  (V-19),  we  see  that  the  quantities  (g  , 

A  2 

(4  me,  .  .  . ,  4  me)  satisfy  the  system  of  equations 


(VI-9) 


■'*cs>  and 


d^cs  k+3 
dt  ”  ^cs 


=  (yMs) 


/  Mc\ 

o  l1  +  mJ 


M  3 
c  res 


__m  _J_ 
M  3 
c  rms 


Me  1  Mm  1  Vk  AVj.  _1  k 

M  3  *  3  Pcs+M,  Mi  3  3  /m< 

c  res  c  rms/  c  c  'res  W 


k  =  1,  2,  3  (VI-10) 
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M  [  e  1  m' 
c 


4  =  4  =  0  when  t  =  t 
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^  me  _  k+3 

dt  “  Mne 


.  M  .  M/_l _ 1_\  k  _  |me 

^  s^o  M  (  3  3  ),vme  3 

rine'  rme 
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„  k+3  ,  . 

=  £  =0  wnen  t  =  t 

^  me  o 


where  pcg  =  J(y^)2  +  (y^)2  +  (y2/  and  pmfi  =  J(y^/  +  (y  2  /  +  (y^/ 


^cs  “  Kiv cs'  T 

6  .  ,  ,  1 
•.V  )  and  (y 


The  quantities 


(y  .  ; . . ,  y  )  and  (y  .  . . . ,  y  )  are  known  as  functions  of  time  from  the  formulas  in  Sec.  II-B 

so  that,  if  we  numerically  integrate  (VI-10)  and  (VI-11)  to  find  (|  ,  . .  . ,  £°s)  and  (l^e-  •  •  •  *  sme) 

the  position  and  velocity  of  the  Earth-Moon  barycenter  and  of  the  Moon  (x  1  ,  . . . ,  x  M  and 
^  ^  CS  cs 

(x  1  ,,,.,x°  )  can  be  determined  from  (VI-9). 

Let  9x^s/9/^  and  9xme/9^in  &  ^  =  1,  ...»  6)  denote  the  partial  derivatives  of  the  position 

and  velocity  of  the  Earth-Moon  barycenter  and  of  the  Moon  with  respect  to  the  initial  osculating 

orbital  elements  (p  \  p  and  (p  1 ,  . . . ,  p  ^  ).  These  quantities  satisfy  differential  equations 
c  c  m  m  .  < 

(V-20)  and  (V-21)  with  initial  conditions  8x^/8 p/  and  8xk  /8 /?•*  (j,  k  =  1,  .  .  . ,  6).  Let 

•  ■  |  UvS  C  Ulllc  XII 

8y  /dp  •*  and  9y  /dp  ^  (j,  k  =  1,  . . . ,  6)  be  the  solutions  of  differential  equations  (11-64)  with 
cs  c  me  in 

these  same  initial  conditions  and  with  p.  =  (yM  )  [1  +  (M  /M  )]  and  p  =  (yM  )  (M  /M  ),  re- 

so  cs  socs 

spectively.  We  define 


j,  k  =  1,  .  . .,  6 


(VI-12) 
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Subtracting  (11-64)  from  (V-20)  and  (V-21),  we  see  that 
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when  t  =  t 


k  =  1,  2,  3 

j=l . 6  .  (VI-14) 


The  quantities  By^Jdpl  and  3yk  /dp}  are  known  as  functions  of  time  from  the  formulas  in 
c  s  c  m  e  m  <  . 

Sec.  II-D  so  that,  if  we  numerically  integrate  (VI-13)  and  (VI-14)  to  find  the  rj  .  and  rj  we  can 

onj 

find  the  partial  derivatives  of  position  and  velocity  with  respect  to  initial  osculating  orbital  ele¬ 
ments  for  the  Earth-Moon  barycenter  and  for  the  Moon  from  formulas  (VI-1?.). 

The  method  of  changing  Encke  orbits  for  the  Earth-Moon  barycenter  and  for  the  Moon  integra¬ 
tions  is  the  same  as  discussed  at  the  end  of  Sec.  VI-A.  It  will  be  necessary  to  change  Encke  orbits 
more  often  in  the  case  of  the  Moon  than  in  the  case  of  the  Earth-Moon  barycenter. 
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APPENDIX  A 

PRECESSION-NUTATION  OF  THE  EARTH 


The  notation  1950.0  denotes  the  beginning  of  the  Besselian  year,  which  is  the  instant  near 
the  beginning  of  the  calendar  year  1950  when  the  right  ascension  of  the  mean  Sun  was  exactly 
18^40m.  In  more  conventional  notation,  1950.0  is  thus  J.E.D.  2433282.423  or  1950  January 

0  .923  Ephemeris  Time  (see  Ref.  5). 

12  3  1 

Let  (X  ,  X  ,  X  )  be  a  rectangular  coordinate  system  whose  X  -axis  points  toward  the  mean 

3 

vernal  equinox  of  1950.0,  whose  X  -axis  points  toward  the  mean  north  pole  of  1950.0,  and  whose 

2  12  3 

X  -axis  completes  the  right-hand  system.  In  more  concise  language,  we  say  that  (X  ,  X  ,  X  ) 

is  a  rectangular  coordinate  system  referred  to  the  mean  equinox  and  equator  of  1950.0  This  is 

12  3 

the  coordinate  system  in  which  we  are  going  to  integrate  the  equations  of  motion.  Let  (x  ,  x  ,  x  ) 

be  a  coordinate  system  referred  to  the  true  equinox  and  equator  of  date  with  the  same  origin  as 
12  3  1 

the  (X  ,  X  ,  X  )  coordinate  system.  The  x  -axis  points  toward  the  true  vernal  equinox  of  date, 

3  2 

the  x  -axis  points  toward  the  true  north  poie  of  date,  and  the  x  -axis  completes  the  right-hand 

12  3  12  3 

system.  The  relation  between  the  (x  ,  x  ,  x  )  and  (X  ,  X  ,  X  )  coordinate  systems  is 


(A-l) 


(A-2) 


with  N  and  P  being  the  nutation  and  precession  matrices,  respectively.  The  matrix  A  appears 
in  formulas  (III-55)  and  (III-79). 

We  now  give  the  established  expressions  for  the  precession  and  nutation.  First,  we  follow 
Ref.  37  in  defining  the  angles 

f  =  2304V948T  +  0'.'302T2  +  0V0179T3 


7,  =  2304V948T  +  l'.'093T2  +  0VO192T3 


0  =  2004'! 255T  -  0'.'426T2  -  0V0416T3  (A-3) 

where  T  is  measured  in  tropical  centuries  of  36524.21988  ephemeris  days  from  the  epoch  1950.0 
(J.E.D.  2433282.423)  to  the  instant  of  interest.  Then  the  precession  matrix  at  this  instant  is 
given  by38 
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P..=cos£  cos  0  cosz  -  sin  £  sinz 

11  o  o 

P.0=-sin£  cos  0  cosz  -  cos  £  sinz 

12  “o  o 


P ,  =  —  sin  0  cos  z 
13 


P21  =  cos  £q  cos  0  sin  z  +  sin  £q  cos  z 
P,9  =  -  sin  £  cos  0  sin  z  +  cos  £  cos  z 


P23  =  -  sin 0  sinz 


P,  .  =  cos  £  sin  © 
31  o 

P^2  =  —  sin  £q  sin  0 


P33  =  cos  0 


(A-4) 


The  mean  obliquity  of  the  ecliptic  is ' 


eQ  =  23°27,08'.,26  -  46'.'845T  -  0'.'0059T2  +  O'.’OOISIT3 


(A-5) 


where  T  is  measured  in  Julian  centuries  of  36525  ephemeris  days  from  the  epoch  1900  January 
0.5  E.T.  =  J.E.D.  2415020.0  to  theinstant  of  interest.  Let  and  Af  be  the  nutations  in  longitude 
and  obliquity,  respectively,  as  given  by  the  series  in  Ref.  40.  The  true  obliquity  of  the  ecliptic 
is  then 


e  =  e  +  A<r 
o 


(A-6) 


Finally,  the  nutation  matrix  is  given  by 


N11  NaZ  N13 


A^cose  -Arsine 


N  =  N2i  N22  N23 


A  ip  cos  e 


—  Ae  .  (A-7) 


N  N  N 

in3i  in32  in33 


A  ip  sin*- 


appendix  b 

ROTATION  AND  PHYSICAL  LIBRATION  OF  THE  MOON 


12  3 

Let  {y  ,  y  ,  y  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Moon  whose 

3 

coordinate  axes  point  along  the  principal  axes  of  inertia  of  the  Moon.  The  y  -axis  points  along 

1 

the  axis  of  rotation  of  the  Moon,  while  the  y  -axis  always  points  in  the  general  direction  of  the 

Earth,  the  period  of  rotation  of  the  Moon  about  its  center  of  mass  being  the  same  as  its  orbital 

2  12 
period..  The  y  -axis  completes  the  right-hand  system,  so  that  the  (y  ,  y  )  plane  is  the  plane  of 

the  Moon's  equator.  In  this  coordinate  system,  the  second  harmonic  of  the  gravitational  potential 

of  the  Moon  has  the  form  (III-70). 

12  3 

Let  (u  ,  u  ,  u  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Moor,  re¬ 
ferred  to  the  mean  equinox  and  ecliptic  of  date,  and  suppose  that  tP  is  the  longitude  of  the  de¬ 
scending  node  of  the  lunar  equator  on  the  ecliptic  of  date  measured  from  the  mean  equinox  of 

date,  that  0  is  the  inclination  of  the  lunar  equator  on  the  ecliptic  of  date,  and  that  ip  is  the  angu- 

1  12  3 

lar  distance  of  the  positive  part  of  the  y  -axis  of  the  coordinate  system  (y  ,  y  ,  y  )  from  the 
descending  node  of  the  lunar  equator.  Then  (II-l)  and  (II-2),  with  fi  =  tp,  i  =  -  0  and  w  =  <p,  imply 


3^  =  Z  V 


j  =  1,  2,  3 


(B-l) 


where 


J  -  2  v' 

<=  1 


U  ^  ^  =  cos  tp  cos  tp  -  sin  tp  sin  tp  cos  © 
U^2  =  sirup  cos  tp  +  cos  tp  sin<p  cos© 
U13  =  —  sin  (p  sin  0 

U2^  =  —  cos  tp  sin  <p  —  sin  tp  cos  tp  cos  0 
U22  =  —  sin  tp  sin  tp  +  cos  tp  cos  tp  cos  0 


«  =  —  cos  tp  sin  0 


U, .  =  -  sin  tp  sin  © 


U^2  =  cos  tp  sin  0 


U33  =  cos  © 


(B-2) 


12  3 

Let  (x  ,  x  ,  x  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Moon  re¬ 
ferred  to  the  mean  equinox  and  equator  (of  the  Earth)  of  date.  Let  e  be  the  mean  inclination 

42  ° 

of  the  ecliptic  as  given  in  (A-5).  Then  we  can  write 
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uj=  Z  VMx* 


j  =  1.  2,  3 


(B-3) 


where 


x3  =  V  V .  .u* 

L  fj 

i  =  1  J 


- 1 

<J 

V12 

V13_ 

1 

0 

0 

V21 

V22 

V23 

= 

0 

cos  e 

o 

sin  e 

o 

731 

V32 

V33_ 

0 

-sin  e 

0 

cos  e 

0 

(B— 4) 


Combining  (B-l)  and  (B-3),  we  see  that 


5^=  I  W„ 


j  =  1,  2,  3 


(B-5) 


where 


*j  =  Z  wf/ 


w.f  =  y  u..  v. . 

jt  L  jk  kf 
k=l 


j.  <  =  1,  2,  3 


(B-6) 


so  that,  by  (B-4),  we  have 

W.,  =  U., 
J1  J1 


W._  =  U.,  cos  e  —  U.-  sin  e  j  =  1,  2,  3 
j2  j2  o  j3  o  J 

W._  =  U._  sine  +  U.,  cose 
j3  j2  o  ]3  o) 


(B-7) 


12  3 

Let  (X  ,  X  ,  X  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Moon  re¬ 
ferred  to  the  mean  equinox  and  equator  (of  the  Earth)  of  1950.0,  the  reference  system  in  which 
we  are  going  to  integrate  the  equations  of  motion.  Then  we  have 


I  P„x' 


j  =  i,  2,  3 


(B-8) 


xj=  Z  P, 
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where  (P^)  is  the  precession  matrix  of  (A-4).  Combining  (B-6)  and  (B-8),  we  have 


y 5  -  Z  Bj(x' 


f=l 

3 


Xj=  E 


f=l 


j  =  1,  2,  3 


(B-9) 


where 


V  ^  wikpk«  •  i-  z’  3  • 

k=l 


(B-10) 


This  is  the  matrix  B  which  appears  in  (III-75)  and  (III-79). 

Let  (  be  the  mean  longitude  of  the  Moon,  measured  in  the  ecliptic  from  the  mean  equinox 

of  date  to  the  mean  ascending  node  of  the  lunar  orbit  and  then  along  the  orbit.  Let  Q  be  the 

longitude  of  the  mean  ascending  node  of  the  lunar  orbit  on  the  ecliptic  measured  from  the  mean 

equinox  of  date.  Finally,  let  I  be  the  inclination  of  the  mean  lunar  equator  to  the  ecliptic.  Then 

43 

the  angles  ip,  @,  cp  of  formulas  (B-2)  are  * 
ij>  =  ft  +  cr 
0  =  I  +  p 

<p  =  180°  +  ((-  ft)  +  (r  -  a)  (B-ll) 

where  a,  p  and  r  are  the  physical  librations  in  node,  inclination  and  longitude,  respectively. 

We  now  determine  the  quantities  on  the  right  side  of  (B-ll).  First,  the  inclination  of  the 

44  45 

mean  lunar  equator  on  the  ecliptic  is  ’ 

I  =  lo32'20"  =  1?  53889 


=  0.0268587  radian  (B-12) 

Next,  according  to  Ref.  >6,  we  have 

Q  =  259? 183275  -  0?0529539222d 

+  1?  557  X  10'12d2  +  5 f  0  X  10'20d3 
(  -  Q  -  Ilf 250889  +  13f2293504490d 

—  2? 407  X  10"12d2  -  If  1  X  10-20d3  (B-13) 

where  d  is  the  number  of  days  that  have  elapsed  from  J.E.D.  2415020.0.  Finally,  the  physical 

47 

libration  of  the  Moon  is 

r  =  —  1 2 V 9  sinf  —  0'!3  sin2f  +  6 5 V 2  si n f 

+  9V 7  sin (2F  -  21)  +  IV 4  sin(2F  -  2D)  +  2V5  sin(D-  f) 

-0V6  sin  (2D  -  U  +  (')  -  7V3  sin  (2D  -  2f) 

—  3V0  sin  (2D  —  t)  —  0V4  sin  2D  +  7V6  sinfl  ;  (B-14) 
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1 


106"  cos  t  4  35"  cos  (2F  —  f)  -  11"  cos  2F 


P  =  - 

-  3"  cos(2F  -  2D)  -  2"  cos  (2D  -  ()  ;  (B-15) 

I(t  -a)  =  108"  sinf-  35"  sin(2F-  f)  4  11"  sin2F 

4  3"  sin  (2F  -  2D)  4  2"  sin  (2D  -  f)  (B-16) 

where  I  is  given  by  (B-12)  measured  in  radians,  where  we  have  taken  the  parameter  f  in  Ref.  47 

to  be  f  =  0.73,  and  where  the  arguments  t,  £',  F  and  D  are  given  in  Ref.  40  as  functions  of  time. 

The  relations  between  the  arguments  £,  £',  F  and  D,  and  the  arguments  g,  g',  w  and  u>'  in  Ref.  47 

.48 

are  given  by 

(  =  g  g  =  l 

('  -  g1  g'  =  (' 

•  .  ( B  - 1 7 ) 

D  =  g1  -  g1  +  u  -  u>’  o)  =  F  —  ( 

F  =  g  4  cj  «'  =  F  —  D  —  f. 
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APPENDIX  C 

ORIENTATION  OF  THE  SUN 


i 


According  to  Ref.  49,  we  have 

Inclination  of  solar  equator  to  ecliptic 

Longitude  of  ascending  node  of  solar 
equator  on  ecliptic 

where  t  is  the  time  in  years  from  1850.  Thus,  in  1950.0,  the  longitude  of  the  ascending  node 
of  the  solar  equator  on  the  ecliptic  was 

SI  =  75°3!75  =  75?0625  .  (C-2) 

s 

Now  the  rate  of  the  precession  of  the  equinox  backward  along  the  ecliptic  is  50".25  per  year.  Since 

formulas  (C  —  1 )  were  derived  from  observations,  we  can  conclude  that  no  precession  of  the  solar 

equator  along  the  ecliptic  has  been  observed.  If  the  Sun  had  an  equatorial  bulge,  such  a  precession 

would  arise  from  the  gravitational  action  of  the  planets.  Thus,  the  fact  that  no  precession  of 

the  solar  equator  has  been  observed  puts  an  upper  limit  on  the  possible  magnitude  of  the  second 

harmonic  of  the  Sun's  gravitational  potential.  However,  the  planetary  torques  acting  on  such  a 

solar  equatorial  bulge  would  be  so  small  that  this  upper  limit  is  not  much  of  a  restriction. 

12  3 

Let  (x  ,  x  ,  x  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Sun  whose 

3  1 

x  -axis  points  toward  the  north  pole  of  the  Sun,  whose  x  -axis  is  the  intersection  of  the  equator 

2 

of  the  Sun  and  the  mean  ecliptic  of  1950.0,  and  whose  x  -axis  completes  the  right-hand  system. 
12  3 

Let  (u  ,  u  ,  u  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Sun  referred 
to  the  mean  equinox  and  ecliptic  of  1950.0.  Then  the  results  of  Sec.  II-A  imply  that 


I  =  7° 15' 
s 


SI  (t)  =  73°40'  +  50'.'25t 
s 


(C-l) 


11  2 
x  =  u  cos  11  +  u  sin  SI 
s  s 

2  1  2  3 

x  = -u  sinS2  cos  I  +u  cos  SI  cos  I  +u  sin  I 
s  s  s  s  s 

3  1  2  3 

x  =  u  sin  SI  sin  I  -u  cos  SI  sin  I  +  u  cos  I 
s  s  s  s  s 


(C-3) 


12  3 

Let  (X  ,  X  ,  X  )  be  the  coordinate  system  with  origin  at  the  center  of  mass  of  the  Sun  referred 

1  2'  3  12  3 

to  the  mean  equinox  and  equator  of  1950.0.  The  relation  between  the  (u  ,  u  '  u  )  and  (X  ,  X  ,  X  ) 

coordinate  systems  is  given  by  (B-3),  with  eQ  =  e  being  the  mean  inclination  of  the  ecliptic  in 

1950.0.  Combining  (C-3)  and  (B-3),  we  see  that 

11  2  _  3  _ 

x  =X  cos  SI  rX  sin  SI  cos  <r  +  X  sin  SI  sine 
s  s  s 

2  12  _ 
x  =  — X  sin  SI  cos  I  +X(cosS2  cos  I  cos  e  —  sin  I  sin  e) 
s  s  s  s  s 

3  _ 

+  X  (cos  SI  cos  I  sin  e  +  sin  I  cos  e) 
s  s  s 

3  1  2  -  - 

x  =X  sin  SI  sin  I  -X  (cos  S2  sin  I  cos  e  +  cos  I  sin  e) 
s  s  s  s  s 


+  X  (—cos  SI  sin  I  sin  e  +  cos  I  cos  e) 
s  s  s 


(C-4) 
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Finally,  comparing  (III-47)  and  (C-4),  we  see  that  the  quantities  C3[.  (k  =  1,  2,  3)  appearing 
in  (III-50)  are 

C31  =  sin  Slg  sin  Ig 

C  =  —  cos  ft  sin  I  cos  e  —  cos  I  sine 
s  s  s 

C-j,  =  —  cos  ft  sin  I  sin  e  +  cos  I  cos?  .  (C-5) 
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